A Reevaluation of the Italian PGA Data in a Varying Coefficient Ground-Motion Model with INLA
Nicolas Kuehn1
10 September, 2021
- 1 Introduction
- 2 Getting Started
- 3 Model Description
- 4 Read Data
- 5 Non-Spatial INLA Models
- 6 Spatial INLA Models
- 6.1 Generate Mesh
- 6.2 Make SPDE Priors and Projector Matrices
- 6.3 Spatial Model with only spatially varying constant terms
- 6.4 Spatial Model with only spatially varying event constant terms
- 6.5 Spatial Model with only spatially varying station constant term
- 6.6 Spatial Model with spatially varying constant and geometrical spreading terms
- 6.7 Spatial Model with spatially varying constant, and anelastic attenuation terms
- 6.8 Spatial Model with spatially varying constant and \(V_{S30}\) terms
- 6.9 Spatial Model with spatially varying constant and all terms
- 6.10 Spatial Model without random effect, spatially varying constant, wth lnR,R,VS
- 6.11 Spatial Model with spatially varying constant and cell-specific attenuation
- 6.12 Spatial Model with all spatially varying terms and cell-specific attenuation
- 6.13 Spatial Model with spatially varying constants and geometrical spreading term and cell-specific attenuation
- 6.14 Spatial Model with spatially varying constants and VS30 scaling term and cell-specific attenuation
- 7 MS-GWR
- 8 Model Comparison
- 9 Comparisons of Spatial Coefficients
- 10 Prediction
- 11 Some Remarks
- 12 Some additional plots
- 13 Acknowledgments
- References
1 Introduction
In recent years, there has been a push to relax the ergodic assumption (Anderson and Brune 1999) in the development of empirical ground-motion models (GMMs). Loosely speaking, the ergodic assumption says that one can trade time for space, and that the ground-motion distribution a a particular site is the same as the ground-motion distribution average over many sites. This means that an ergodic GMM ignores repeatable systematic source, path, and site effects. Due to increasing numbers of observed data, it is possible to estimate these effects. Much has been written about nonergodic GMMs and their use in probabilistic seismic hazard analysis (PSHA). For a schematic overview of nonergodic PSHA, see e.g. Villani and Abrahamson (2015) or Walling and Abrahamson (2012). Abrahamson et al. (2019) conducted a nonergodic PSHA for three sites in California, and showed that there can be large differences in PSHA results due to the inclusion of nonergodic effects on ground motion.
Landwehr et al. (2016) introduced the concept of spatially varying coefficient models (VCMs) (Gelfand et al. 2003; Bussas et al. 2017) to the estimation of nonergodic GMMs. In a VCM, several coefficients are functions of location, thus encoding spatial variation of ground-motion scaling. In the model of Landwehr et al. (2016), the spatially varying coefficients are assumed to be spatially correlated, which enables their estimation, as well as extrapolation to new location, More precisely, the varying coefficients are assumed to be distributed according to a Gaussian process (GP) (Rasmussen and Williams 2006). The GP can be understood as a prior distribution over the functional form of the spatial dependence of the coefficients, with the covariance function providing some constraints on the shape of the function. In Landwehr et al. (2016), the model was estimated using the Matlab GPML package (Rasmussen and Nickisch 2010), using the FITC approximation (Snelson and Ghahramani 2005) to the covariance matrix to make the model computationally tractable. Kuehn (2021) showed how VCMs based on GPs can be estimated via Bayesian inference, using Markov Chain Monte Carlo (MCMC) sampling using he computer language Stan (Carpenter et al. 2017), and the integrated nested Laplace approximation INLA) (Rue, Martino, and Chopin 2009). Results showed that both approaches give comparable results, with INLA being faster. Several studies have estimated nonergodic GMMs based on GP VCMs (Lavrentiadis, Abrahamson, and Kuehn 2021) have been estimated recently.
Recently, Caramenti et al. (2020) introduced an alternative to GPs for the estimation of spatially varying GMMs, based on geographically weighted regression (GWR) (Fotheringham, Brunsdon, and Charlton 2002), called multi-source geographically weighted regression (MS-GWR). They showed an example based on Italian strong motion (PGA) data, using the ITA2018 GMM (Lanzano et al. 2019) as a reference. Lanzano et al. (2021) investigated MS-GWR further, extending the model to different response spectral periods. They found that the MS-GWR model provides an improvement in terms of predictive accuracy over a linear base model
GWR and GP based VCMs are based on different underlying assumptions, which naturally raises the question how they compare. In GMM development, epistemic uncertainty is large (in particular for nonergodic models), so there is value in having different models based on different methods. GWR constructs local models for different locations, based on weighting data according to their distance to the location. The weights are determined by a spatial kernel. By contrast, GP based models assume an underlying functional dependence of the spatial coefficients, and place a GP prior on the functions.
Here, we fit different spatial models with INLA, and compare the fit to an ergodic base model, as well as the model of Caramenti et al. (2020).
2 Getting Started
First, we load the packages required for the analysis.
3 Model Description
The functional form of the ITA18 model (Lanzano et al. 2019) is \[ \begin{aligned} \log_{10} PGA &= a + b_1 (M_W - M_h) \; {1}_{(M_w \leq M_h)} + b_2 (M_W - M_h) \; {1}_{(M_w > M_h)} \\ &\quad + \left[ c_2 + c_1 (M_W - M_h)\right] \log_{10} \sqrt{R_{JB}^2 + h^2} + c_3 \sqrt{R_{JB}^2 + h^2} \\ &\quad + k \left[ \log_{10} \frac{V_{S30}}{800} \; {1}_{(V_{S30} \leq 1500)} + \frac{1500}{800} \; {1}_{(V_{S30} > 1500)}\right] + f_1 F_{SS} + f_2 F_R \\ &\quad + \delta B + \delta S2S + \delta WS \end{aligned} \] where \(1_{(M_w \leq M_h)}\) is an indicator function that evaluates to one if the condition is true, and zero otherwise. \(F_{SS}\) and \(F_R\) are indicators for strike-slip and reverse-faulting event, respectively. We use the same functional form, and fit different models to the PGA data used in Caramenti et al. (2020). For the spatial models, we make the constant \(a\), the geometrical spreading coefficient \(c_2\), the anelastic attenuation coefficient \(c_3\), and the \(V_{S30}\)-scaling coefficient \(k\) spatially varying. There are two spatially varying constants, one dependent on the event coordinate \(\vec{x}_e\), the other dependent on the station coordinate \(\vec{x}_s\). The spatially varying geometrical spreading coefficient depends on \(\vec{x}_e\), the \(V_{S30}\)-coefficient on \(\vec{x}_s\). Following Landwehr et al. (2016), Caramenti et al. (2020) modeled the anelastic attenuation coefficient as dependent on the source coordinate \(\vec{x}_e\). We do the same here, but additionally test a model where the anelastic attenuation coefficient depends n the midpoint between the event and station coordinate, \(\vec{x}_{es}\). All spatial terms are modeled as adjustment terms to the fixed effects. Thus, the full spatial model becomes \[ \begin{aligned} \log_{10} PGA &= a + \delta a_e(\vec{x}_e) + \delta a_s(\vec{x}_s) + b_1 (M_W - M_h) \; {1}_{(M_w \leq M_h)} + b_2 (M_W - M_h) \; {1}_{(M_w > M_h)} \\ &\quad + \left[ (c_2 + \delta c_2(\vec{x}_e)) + c_1 (M_W - M_h)\right] \log_{10} \sqrt{R_{JB}^2 + h^2} + (c_3 + \delta c_3(\vec{x}_{e/es})) \sqrt{R_{JB}^2 + h^2} \\ &\quad + (k + \delta k(\vec{x}_s)) \left[ \log_{10} \frac{V_{S30}}{800} \; {1}_{(V_{S30} \leq 1500)} + \frac{1500}{800} \; {1}_{(V_{S30} > 1500)}\right] + f_1 F_{SS} + f_2 F_R \\ &\quad + \delta B_0 + \delta S2S_0 + \delta WS_0 \end{aligned} \] where we have added the subscript \(_0\) to the event and station term to indicate that they are the terms after the systematic, repeatable nonergodic terms are taken out. We test different combinations of the spatial terms, to see which term gives the best improvement. We also model the anelastic attenuation term using the cell-specific attenuation model (Dawood and Rodriguez-Marek 2013; Kuehn, Abrahamson, and Walling 2019). Compared to a varying coefficient \(c_3\), which only depends on a point coordinate, this has the advantage that it models the anelastic attenuation along a path, which makes more physical sense. In that model, the anelastic attenuation model becomes \[ \begin{aligned} f_{attn}(R) &= c_3 \sqrt{R_{JB}^2 + h^2} + \sum_{i = 1}^{N_C} \delta c_{3,i} \; \Delta R_i \\ & = c_3 \sqrt{R_{JB}^2 + h^2} + \overrightarrow{\Delta R} \cdot \vec{\delta c_3} \end{aligned} \] where \(\delta c_{3,i}\) is the cell-specific adjustment coefficient for the \(i\)th cell, and \(\Delta R_i\) is the fraction of path in the \(i\)th cell. The distances within each cell are calculated based on epicentral distance, since the closest points on the rupture are not available.
The adjustment coefficients are modeled as a Gaussian process (or a Gaussian random field) with mean zero and Mat'ern covariance function \[ \begin{aligned} \delta a &= GP(0, \omega_a^2 k(\vec{x}, \vec{x}')) \\ k(\vec{x}, \vec{x}') &= \frac{2^{(1 - \nu)}}{\Gamma(\nu)} (\kappa |\vec{x} - \vec{x}'|)^\nu K_\nu (\kappa |\vec{x} - \vec{x}'|) \end{aligned} \] where \(\Gamma\) is the gamma function, \(K_\nu\) is the modified Bessel function of the second kind, \(\kappa\) is a scale parameter and \(\nu\) is a smoothness parameter. The cell-specific attenuation coefficients are modeled as random effects, with a normal distribution with mean zero and standard deviation \(\omega_{c3}\) as the prior (following Kuehn, Abrahamson, and Walling (2019)) \[ \delta c_3 \sim N(0, \omega_{c3}) \]
All models are estimated with the R-INLA package (https://www.r-inla.org/) (Rue et al. 2017). INLA provides a computationally efficient way to estimate spatial models based on Gaussian random fields using the stochastic partial differential equations (SPDE) approach (Lindgren, Rue, and Lindström 2011; Lindgren and Rue 2015; Bakka et al. 2018). For \(\nu = 1\) in the Mat'ern covariance function, the Gaussian field can be represented as a Gaussian Markov field (GMRF) by representing a solution to the SPDE using finite elements. This leads to a sparse representation of the precision matrix, which allows for efficient computation. The Gaussian field is approximated using basis functions, which are evaluated on nodes of a triangular mesh (Lindgren, Rue, and Lindström 2011; Krainski et al. 2019).
4 Read Data
Below, we read in the data. The data is taken from the electronic supplement of Caramenti et al. (2020), available at https://github.com/lucaramenti/ms-gwr/. Most of the code is taken from there, reading in the shape file for Italy, and converting coordinates to UTM, as well as preparing the linear covariates. In total, there are 4784 records from 137 events and 927 stations.
dir_esupp <- '/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/Caramenti_esupp/'
dataset = readRDS(file.path(dir_esupp,"italian_data_pga.RData"))
shape_utm = st_read(file.path(dir_esupp,'confini_ut33.shp'))## Reading layer `confini_ut33' from data source
## `/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/Caramenti_esupp/confini_ut33.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -158718.7 ymin: 3930217 xmax: 800075.7 ymax: 5219220
## Projected CRS: WGS 84 / UTM zone 33N#
shape_utm_no_lamp = shape_utm
shape_utm_no_lamp$geometry[[1]][[12]] = NULL
shape_utm_no_sardinia = shape_utm
shape_utm_no_sardinia$geometry[[1]][[3]] = NULL
for (i in 4:10){
shape_utm_no_sardinia$geometry[[1]][[4]] = NULL
}
for (i in 12:57){
shape_utm_no_sardinia$geometry[[1]][[5]] = NULL
}
#Project coordinates of the events to work with UTM33:
latitude_ev = dataset$ev_latitude
longitude_ev = dataset$ev_longitude
long_lat_ev = cbind(longitude_ev, latitude_ev)
utm_ev = project(long_lat_ev, "+proj=utm +zone=33 ellps=WGS84")
utm_ev = as.data.frame(utm_ev)
long_lat_ev = as.data.frame(long_lat_ev)
#Project coordinates of the stations to work with UTM33:
latitude_st = dataset$st_latitude
longitude_st = dataset$st_longitude
long_lat_st = cbind(longitude_st, latitude_st)
utm_st = project(long_lat_st, "+proj=utm +zone=33 ellps=WGS84")
utm_st = as.data.frame(utm_st)
long_lat_st = as.data.frame(long_lat_st)
#Build spatial points data frame for UTM33 coordinates:
utm_ev_sp = SpatialPointsDataFrame(utm_ev, dataset[,1:6])
utm_st_sp = SpatialPointsDataFrame(utm_st, dataset[,1:6])
#Transform shapefile into spatial dataframe:
shape_utm_spatial = as_Spatial(shape_utm_no_sardinia)
grid_utm = makegrid(shape_utm_spatial, cellsize = 10000) # cellsize in map units!
grid_utm = SpatialPoints(grid_utm, proj4string = CRS(proj4string(shape_utm_spatial)))## Warning in proj4string(shape_utm_spatial): CRS object has comment, which is lost
## in outputgrid_inside_utm = grid_utm[shape_utm_spatial,]
coords_utm = grid_inside_utm@coords
coords_df_utm = as.data.frame(coords_utm)
# read n data file with event and station ids
# correcte for Vs30 differences between stations of same location
data <- read.csv(file.path(dir_esupp,'italian_data_pga_id_utm_stat.csv'))
#### Cell-specific distances are computed outside of R
data_dm <- rstan::read_rdump('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/dm_25x25.Rdata')
# Set linear predictors
mh = 5.5
mref = 5.324
h = 6.924
attach(data)
b1 = (mag-mh)*(mag<=mh)
b2 = (mag-mh)*(mag>mh)
c1 = (mag-mref)*log10(sqrt(JB_complete^2+h^2))
c2 = log10(sqrt(JB_complete^2+h^2))
c3 = sqrt(JB_complete^2+h^2)
f1 = as.numeric(fm_type_code == "SS")
f2 = as.numeric(fm_type_code == "TF")
k = log10(vs30/800)*(vs30<=1500)+log10(1500/800)*(vs30>1500)
y = log10(rotD50_pga)
detach(data)
eq <- data$EQID
stat <- data$STATID
n_rec <- length(b1)
n_eq <- max(eq)
n_stat <- max(stat)
n_cell <- data_dm$NCELL
data_reg <- data.frame(Y = y,
intercept = 1,
M1 = b1,
M2 = b2,
MlnR = c1,
lnR = c2,
R = c3,
Fss = f1,
Frv = f2,
lnVS = k,
eq = eq,
stat = stat
)
data_reg$idx_cell <- 1:nrow(data_reg)
data_reg$idx <- 1:n_recggplot() +
geom_point(dataset, mapping = aes(x = JB_complete, y = mag)) +
scale_x_continuous(trans='log10') +
labs(x = "R_JB (km)", y = "M") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
ggplot() +
geom_sf(data = shape_utm_no_lamp, size = 1.6, color = "black", fill = NA )+
geom_point(data = utm_st, aes(x=longitude_st, y=latitude_st), fill= 'firebrick3',
size = 2, shape = 24, stroke = 0.7)+
geom_point(data = utm_ev, aes(x=longitude_ev, y=latitude_ev), fill= 'cadetblue2',
size = 2, shape = 21, stroke = 1.5)+
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Scatterplot of Joyner-Boore distance and magnitude, and map of events and stations.
ggplot() +
geom_point(dataset, mapping = aes(x = mag, y = rotD50_pga, colour = JB_complete)) +
scale_y_continuous(trans='log10') +
labs(x = "M", y = "PGA (cm/s)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
ggplot() +
geom_point(dataset, mapping = aes(x = JB_complete, y = rotD50_pga, colour = mag)) +
scale_y_continuous(trans='log10') + scale_x_continuous(trans='log10') +
labs(x = "R_JB (km)", y = "PGA (cm/s)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)## Warning: Transformation introduced infinite values in continuous x-axis
Scatterplot of PGA vs magnitude and Joyner-Boore distance.
5 Non-Spatial INLA Models
First, we fit an ergodic base model as reference. This model corresponds to the ITA18 model of Lanzano et al. (2019), and is a fit to the base functional form without any spatial terms. This model includes random effects for events and stations, i.e. event terms \(\delta B\) and station terms \(\delta S\). We also fit a model without random effects. In its current form, MS-GWR does not include random effects, so we also have this model as reference.
# priors for standard deviation parameters
prior_prec_tau_lg <- list(prec = list(prior = "loggamma", param = c(0.9, 0.007)))
prior_prec_phiS2S_lg <- list(prec = list(prior = "loggamma", param = c(0.9, 0.007)))
prior_prec_phiSS_lg <- list(prec = list(prior = "loggamma", param = c(0.9, 0.007)))
# priors for fixed effects
prior.fixed <- list(mean.intercept = 0, prec.intercept = 0.04,
mean = (list(R=-0.005, Fss = 0.02, Frv = 0.1, default=0)),
prec = (list(R=10000, Fss = 1500, Frv = 625, default = 0.04)))
form <- Y ~ M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg)
fit_inla <- inla(form,
data = data_reg,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE)
)
# model with only event random effects
form_eq <- Y ~ M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg)
fit_inla_eq <- inla(form_eq,
data = data_reg,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE)
)
# also fit a model without random effects,
# since Caramenti et al (2020) does not include random effects
fit_inla_noranef <- inla(Y ~ M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS,
data = data_reg,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE)
)
df_coeff <- data.frame(mean_M1 = fit_inla$summary.fixed$mean, sd_M1 = fit_inla$summary.fixed$sd,
mean_M2 = fit_inla_eq$summary.fixed$mean, sd_M2 = fit_inla_eq$summary.fixed$sd,
mean_M3 = fit_inla_noranef$summary.fixed$mean, sd_M3 = fit_inla_noranef$summary.fixed$sd,
mean_ITA18 = c(3.4210, 0.1940, -0.0220, -1.4056, 0.2871, -0.0029, 0.0860, 0.0105, -0.3946))
row.names(df_coeff) <- c("a", "b1","b2","c2", "c1", "c3", "f1","f2","k")
knitr::kable(df_coeff, digits = 3, row.names = TRUE)| mean_M1 | sd_M1 | mean_M2 | sd_M2 | mean_M3 | sd_M3 | mean_ITA18 | |
|---|---|---|---|---|---|---|---|
| a | 3.407 | 0.049 | 3.609 | 0.055 | 3.549 | 0.047 | 3.421 |
| b1 | 0.194 | 0.039 | 0.253 | 0.049 | 0.294 | 0.034 | 0.194 |
| b2 | 0.006 | 0.073 | 0.091 | 0.091 | 0.131 | 0.038 | -0.022 |
| c2 | -1.399 | 0.030 | -1.533 | 0.033 | -1.483 | 0.035 | -1.406 |
| c1 | 0.287 | 0.014 | 0.249 | 0.018 | 0.219 | 0.019 | 0.287 |
| c3 | -0.003 | 0.000 | -0.002 | 0.000 | -0.003 | 0.000 | -0.003 |
| f1 | 0.058 | 0.021 | 0.041 | 0.022 | 0.048 | 0.013 | 0.086 |
| f2 | 0.028 | 0.026 | 0.050 | 0.028 | -0.016 | 0.011 | 0.011 |
| k | -0.420 | 0.045 | -0.365 | 0.023 | -0.402 | 0.025 | -0.395 |
6 Spatial INLA Models
For the spatial INLA models, we divide the UTM coordinates by 1000, so that we can wok on a scale of kilometers. This is just a matter of convenience. In INLA, the spatial field s represented as basis functions on a triangular mesh, and projected onto the data locations via a projector matrix \(\textbf{A}\). We need to define the mesh and the projector matrices. We also need to define prior distributions for the spatial hyperparameters. These are based on penalized complexity (PC) priors (Simpson et al. 2017; Fuglstad et al. 2019; Franco-Villoria, Ventrucci, and Rue 2018).
Conceptually, a spatial random effects model looks as in the Equation below. The vector of target variables is distributed according to a normal distribution with mean \(\vec{\mu}\) and standard deviation \(\phi_0\). The mean is composed of the base model (the fixed effects) and the random effects, which means that the likelihood factorizes, and all observations are independent (on could also marginalize over the random effects, i.e. put them into the covariance matrix).
\[ \begin{aligned} \vec{Y} &\sim N(\vec{\mu}, \phi_0) \\ \vec{\mu} =& \vec{\mu}_{base} + \vec{\delta}_{adj} \\ \vec{\delta}_{adj} &\sim GP(0, k(\cdot,\cdot)) \end{aligned} \] With the SPDE approach on a mesh, the model is as below. Now, the random effects are defined on the nodes of the mesh, and are distributed according to a multivariate normal distribution with mean zero and precision matrix \(\textbf{P}_{mesh}\), which has size \(N_{mesh} \times N_{mesh}\), where \(N_{mesh}\) is the number of nodes in the mesh. Since the field is Markov, \(\textbf{P}_{mesh}\) is sparse, and inference on the random effects can be done efficiently. The matrix \(\textbf{A}\) is used to project the random effects to the observation locations.
\[ \begin{aligned} \vec{Y} &\sim N(\vec{\mu}, \phi_0) \\ \vec{\mu} =& \vec{\mu}_{base} + \textbf{A} \; \vec{\delta}_{adj,mesh} \\ \vec{\delta}_{adj,mesh} &\sim GP(0, \textbf{P}^{-1}_{mesh}) \end{aligned} \]
6.1 Generate Mesh
First we generate the triangular mesh. We use the same mesh for event-specific and station-specific terms, and thus combine the coordinates for mesh generation. We use a maximum triangle edge length of 10km, and an outer boundary of 50km.
co_eq <- utm_ev/1000
co_stat <- utm_st/1000
co_eq_stat <- cbind(rowMeans(cbind(co_eq$longitude_ev, co_stat$longitude_st)),
rowMeans(cbind(co_eq$latitude_ev, co_stat$latitude_st)))
coords <- unique(rbind(as.matrix(co_eq), as.matrix(co_stat)))
max.edge2 <- 10
bound.outer2 <- 50
mesh = inla.mesh.2d(loc=as.matrix(coords),
max.edge = c(1,5)*max.edge2,
# - use 5 times max.edge in the outer extension/offset/boundary
cutoff = max.edge2,
offset = c(5 * max.edge2, bound.outer2))
print(paste0("numbr of mesh nodes: ", mesh$n))## [1] "numbr of mesh nodes: 7505"ggplot() + theme_bw() + gg(mesh) +
geom_point(data = as.data.frame(co_stat), aes(x=longitude_st,y=latitude_st), color = set1[1]) +
#geom_point(data = as.data.frame(co_eq_stat), aes(x=V1,y=V2), color = set1[1]) +
geom_point(data = as.data.frame(co_eq), aes(x=longitude_ev,y=latitude_ev), color = set1[2]) +
labs(x="X (km)", y="Y (km)") +
theme(axis.title = element_text(size=30), axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20))
Mesh used for spatial INLA models.
6.2 Make SPDE Priors and Projector Matrices
Now we define the prior distributions for the spatial hyperparameters. As stated before, they are based on PC priors. The PC prior for a 2D-Mat`ern field is set by specifying tail probabilities for the spatial range \(\ell\) and the standard deviation of the field \(\omega\). The spatial range is defined as \(\ell = \frac{\sqrt{8 \nu}}{\kappa}\), and corresponds to the distance where the correlation takes a value of 0.140. The PC prior requires one to specify values \(\ell_0, p_\ell\) and \(\omega_0, p_\omega\) such that \[ \begin{aligned} P(\ell < \ell_0) = p_\ell \\ P(\omega > \omega_0) = p_\omega \end{aligned} \] i.e. we need to specify our prior belief that the range is smaller than a value \(\ell_0\), and the standard deviation is larger than a value \(\omega_0\). We can specify different prior distributions for the different spatial terms. Since we do not have strong prior information, we use the same prior for the range for all terms, and only adjust the prior for the standard deviation for the anelastic attenuation term, and the prior becomes \[ \begin{aligned} P(\ell < 30) = 0.5 \\ P(\omega > 0.3) = 0.01 \\ P(\omega > 0.1) = 0.01 \end{aligned} \] which means we assign a prior probability of 50% that the range is smaller than 30km, and a prior probability of 1% that the standard deviation of the systematic spatial terms is smaller than 0.3 (or 0.1 for the anelastic attenuation coefficient).
spde_stat <- inla.spde2.pcmatern(
# Mesh and smoothness parameter
mesh = mesh, alpha = 2,
# P(practic.range < 0.3) = 0.5
prior.range = c(30, 0.5),
# P(sigma > .3) = 0.01
prior.sigma = c(.3, 0.01))
spde_eq <- inla.spde2.pcmatern(
# Mesh and smoothness parameter
mesh = mesh, alpha = 2,
# P(practic.range < 30) = 0.5
prior.range = c(30, 0.5),
# P(sigma > .3) = 0.01
prior.sigma = c(.3, 0.01))
spde_vs <- inla.spde2.pcmatern(
# Mesh and smoothness parameter
mesh = mesh, alpha = 2,
# P(practic.range < 30) = 0.5
prior.range = c(30, 0.5),
# P(sigma > .3) = 0.01
prior.sigma = c(.3, 0.01))
spde_lnR <- inla.spde2.pcmatern(
# Mesh and smoothness parameter
mesh = mesh, alpha = 2,
# P(practic.range < 30) = 0.5
prior.range = c(30, 0.5),
# P(sigma > .3) = 0.01
prior.sigma = c(.3, 0.01))
spde_R <- inla.spde2.pcmatern(
# Mesh and smoothness parameter
mesh = mesh, alpha = 2,
# P(practic.range < 30) = 0.5
prior.range = c(30, 0.5),
# P(sigma > .1) = 0.01
prior.sigma = c(.01, 0.01)) Next, we specify the projector matrices \(\textbf{A}\), which connect the mesh nodes with the data locations. The data coordinates are either the event or station locations (or the midpoint for the anelastic attenuation). For the spatial coefficients associated with covariates, we need to specify these as well. We also create the data stack, which combines the information for all data and effects that go into the model.
A_stat <- inla.spde.make.A(mesh, loc = as.matrix(co_stat))
A_eq <- inla.spde.make.A(mesh, loc = as.matrix(co_eq))
A_lnr <- inla.spde.make.A(mesh, loc = as.matrix(co_eq),
weights = data_reg$lnR)
# linear R term based on event coordinates
A_r <- inla.spde.make.A(mesh, loc = as.matrix(co_eq),
weights = data_reg$R)
# linear R term based on midpoint between event and station coordinates
A_r_es <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_stat),
weights = data_reg$R)
A_vs <- inla.spde.make.A(mesh, loc = as.matrix(co_stat),
weights = data_reg$lnVS)
idx.stat <- inla.spde.make.index("idx.stat",spde_stat$n.spde)
idx.eq <- inla.spde.make.index("idx.eq",spde_eq$n.spde)
idx.lnr <- inla.spde.make.index("idx.lnr",spde_lnR$n.spde)
idx.r <- inla.spde.make.index("idx.r",spde_R$n.spde)
idx.vs <- inla.spde.make.index("idx.vs",spde_vs$n.spde)
# create the stack
stk_spatial_full <- inla.stack(
data = list(y = data_reg$Y),
A = list(A_eq, A_lnr, A_r_es, A_stat, A_vs, 1),
effects = list(idx.eq = idx.eq,
idx.lnr = idx.lnr,
idx.r = idx.r,
idx.stat = idx.stat,
idx.vs = idx.vs,
data_reg
),
tag = 'model_spatial_full')6.3 Spatial Model with only spatially varying constant terms
Here, we fit the reference spatial model, with only spatial constants (\(\delta a_e(\vec{x}_e)\) and \(\delta a_e(\vec{x}_e)\)). We define the formula of the model, which includes the fixed effects, iid random effects for the non-spatially event and station terms (\(\delta B\) and \(\delta S2S\)), and the spatially varying terms (using the previously defined SPDE models).
form_spatial <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq)
fit_inla_spatial <- inla(form_spatial,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.4 Spatial Model with only spatially varying event constant terms
form_spatial_eq <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.eq, model = spde_eq)
fit_inla_spatial_eq <- inla(form_spatial_eq,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.5 Spatial Model with only spatially varying station constant term
form_spatial_stat <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat)
fit_inla_spatial_stat <- inla(form_spatial_stat,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.6 Spatial Model with spatially varying constant and geometrical spreading terms
form_spatial_gs <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.lnr, model = spde_lnR)
fit_inla_spatial_gs <- inla(form_spatial_gs,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.7 Spatial Model with spatially varying constant, and anelastic attenuation terms
Here we estimate models with spatially varying event and station constants, and a spatially varying anelastic attenuation coefficient. Two models are estimated, one where the attenuation depends on the midpoint of event/station coordinates, and one where it depends on event location. The stack that we defined before, stk_spatial_full, includes the projector matrix for the midpoint, so we define a separate stack for the second model.
form_spatial_r <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.r, model = spde_R)
fit_inla_spatial_r <- inla(form_spatial_r,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)# create the stack
stk_spatial_r2 <- inla.stack(
data = list(y = data_reg$Y),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_eq, A_stat, A_r, 1),
effects = list(idx.eq = idx.eq,
idx.stat = idx.stat,
idx.r = idx.r,
data_reg
),
tag = 'model_spatial_r2')
fit_inla_spatial_r2 <- inla(form_spatial_r,
data = inla.stack.data(stk_spatial_r2),
control.predictor = list(A = inla.stack.A(stk_spatial_r2)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.8 Spatial Model with spatially varying constant and \(V_{S30}\) terms
form_spatial_vs <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.vs, model = spde_vs)
fit_inla_spatial_vs <- inla(form_spatial_vs,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.9 Spatial Model with spatially varying constant and all terms
form_spatial_full <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.lnr, model = spde_lnR) + f(idx.r, model = spde_R) +
f(idx.vs, model = spde_vs)
fit_inla_spatial_full <- inla(form_spatial_full,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.10 Spatial Model without random effect, spatially varying constant, wth lnR,R,VS
form_spatial_full_b <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(idx.lnr, model = spde_lnR) + f(idx.r, model = spde_R) +
f(idx.vs, model = spde_vs)
fit_inla_spatial_full_b <- inla(form_spatial_full_b,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)6.11 Spatial Model with spatially varying constant and cell-specific attenuation
Now we fit a model with cell-specific attenuation. This is modeled as a generic Z model in the call to inla, as the cell-specific attenuation terms for all records can be written as \(\textbf{R}_C \; \vec{\delta c}_3\), where \(\textbf{R}_C\) is an \(N \times N_C\) matrix whose rows correspond to the records, columns correspond to the cells, and where \(\textbf{R}_{C,ij}\) is the length of the path in cell \(j\) for record \(i\). For the standard deviation \(\omega_{c3}\) of the cell-specific attenuation coefficients, we use again a PC prior. From Kuehn, Abrahamson, and Walling (2019), the standard deviation of the cell-specific attenuation coefficients in California has a value of about \(\omega_{c3, calif} \approx 0.0054\), so we set the prior assuming that a value larger than 0.01 is unlikely
\[ P(\omega_{c3} > 0.01) = 0.01 \]
dm_sparse <- as(data_dm$RC,"dgCMatrix")
prior_prec_cell <- list(prec = list(prior = 'pc.prec', param = c(0.01, 0.01)))
form_spatial_cell <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
fit_inla_spatial_cell <- inla(form_spatial_cell,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'z' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.We also estimate this model without the empirical Bayes approximation.
fit_inla_spatial_cell_int <- inla(form_spatial_cell,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg))
)## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'z' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.6.12 Spatial Model with all spatially varying terms and cell-specific attenuation
form_spatial_full_cell <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.lnr, model = spde_lnR) + f(idx.vs, model = spde_vs) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
# set starting values based on previous estimated models
theta <- c(fit_inla_spatial_full$mode$theta[c(1,2,3,4,5,6,7,8,9,12,13)],fit_inla_spatial_cell$mode$theta[8])
fit_inla_spatial_full_cell <- inla(form_spatial_full_cell,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)6.13 Spatial Model with spatially varying constants and geometrical spreading term and cell-specific attenuation
form_spatial_gs_cell <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.lnr, model = spde_lnR) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
# set starting values based on previous estimated models
theta <- c(fit_inla_spatial_full$mode$theta[c(1,2,3,4,5,6,7,8,9)],fit_inla_spatial_cell$mode$theta[8])
fit_inla_spatial_gs_cell <- inla(form_spatial_gs_cell,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)6.14 Spatial Model with spatially varying constants and VS30 scaling term and cell-specific attenuation
form_spatial_vs_cell <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.vs, model = spde_vs) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
# set starting values based on previous estimated models
theta <- c(fit_inla_spatial_full$mode$theta[c(1,2,3,4,5,6,7,12,13)],fit_inla_spatial_cell$mode$theta[8])
fit_inla_spatial_vs_cell <- inla(form_spatial_vs_cell,
data = inla.stack.data(stk_spatial_full),
control.predictor = list(A = inla.stack.A(stk_spatial_full)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)7 MS-GWR
Now, we fit the MS-GWR model. This is taken directly from the https://github.com/lucaramenti/ms-gwr/. The calculations are done outside of his document, and read in.
# best bandwidth is set
bwe = 25000
bws = 75000
coords_df_utm = as.data.frame(coords_utm)
sec = SEC_only_calibration(Xc, Xe, Xs, y, "c", bwe, bws, coordinates(utm_ev_sp), coordinates(utm_st_sp))
#Now we can compute all the regression coefficients:
result = SEC_grid_creation(Xc, Xe, Xs, y,"c", bwe, bws, coordinates(utm_ev_sp),
coordinates(utm_st_sp), coords_utm, sec)load(file = file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/','RESULTS', 'results_caramenti.Rdata'))
beta_const = result$beta_c
beta_k = t(result$beta_s)
beta_k_coord = cbind(coords_utm, beta_k)
beta_k_coord = as.data.frame(beta_k_coord)
beta_c2 = result$beta_e[1,]
beta_c2_coord = cbind(coords_utm, beta_c2)
beta_c2_coord = as.data.frame(beta_c2_coord)
beta_c3 = result$beta_e[2,]
beta_c3_coord = cbind(coords_utm, beta_c3)
beta_c3_coord = as.data.frame(beta_c3_coord)
time_msgwr <- 14189.787load(file = file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/',
'RESULTS', 'residuals_caramenti.Rdata'))
data_reg$resid <- data_reg$Y - prediction_msgwr$y0
fit_inla_resid <- inla(resid ~ 1 +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg),
data = data_reg,
family="gaussian",
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE)
)8 Model Comparison
Now that we have fit many different models, we need to compare them. We calculated DIC (Spiegelhalter et al. 2002) and WAIC (Watanabe 2013) for the models, which we will use for model comparison. WAIC and DIC are an approximation of the generalization error, and can be (loosely) thought of as Bayesian extensions of the Akaike Information criterion (AIC) (Akaike 1973). One should note that WAIC ad DIC wok at the observation level, i.e. they contain information about the random effects. For application of GMMs, we typically want to evaluate the models on new events. To do that, we also perform a cross-validation (CV), leaving out all observations from events, training the model on the remaining data, and evaluating the models on the left out events.
8.1 Model Comparison based on Information Criteria
Here we compare all the models fit using INLA via their WAIC and DIC values. We also compare the CPU time (not important for model comparison, but important for applicability).
df_comp <- data.frame(matrix(c(fit_inla$waic$waic, fit_inla$dic$dic, fit_inla$cpu.used[4],
fit_inla_eq$waic$waic, fit_inla_eq$dic$dic, fit_inla_eq$cpu.used[4],
fit_inla_noranef$waic$waic, fit_inla_noranef$dic$dic, fit_inla_noranef$cpu.used[4],
fit_inla_spatial$waic$waic, fit_inla_spatial$dic$dic, fit_inla_spatial$cpu.used[4],
fit_inla_spatial_eq$waic$waic, fit_inla_spatial_eq$dic$dic, fit_inla_spatial_eq$cpu.used[4],
fit_inla_spatial_stat$waic$waic, fit_inla_spatial_stat$dic$dic, fit_inla_spatial_stat$cpu.used[4],
fit_inla_spatial_gs$waic$waic, fit_inla_spatial_gs$dic$dic, fit_inla_spatial_gs$cpu.used[4],
fit_inla_spatial_r$waic$waic, fit_inla_spatial_r$dic$dic, fit_inla_spatial_r$cpu.used[4],
fit_inla_spatial_r2$waic$waic, fit_inla_spatial_r2$dic$dic, fit_inla_spatial_r2$cpu.used[4],
fit_inla_spatial_vs$waic$waic, fit_inla_spatial_vs$dic$dic, fit_inla_spatial_vs$cpu.used[4],
fit_inla_spatial_full$waic$waic, fit_inla_spatial_full$dic$dic, fit_inla_spatial_full$cpu.used[4],
fit_inla_spatial_full_b$waic$waic, fit_inla_spatial_full_b$dic$dic, fit_inla_spatial_full_b$cpu.used[4],
fit_inla_spatial_cell$waic$waic, fit_inla_spatial_cell$dic$dic, fit_inla_spatial_cell$cpu.used[4],
fit_inla_spatial_cell_int$waic$waic, fit_inla_spatial_cell_int$dic$dic, fit_inla_spatial_cell_int$cpu.used[4],
fit_inla_spatial_full_cell$waic$waic, fit_inla_spatial_full_cell$dic$dic, fit_inla_spatial_full_cell$cpu.used[4],
fit_inla_spatial_gs_cell$waic$waic, fit_inla_spatial_gs_cell$dic$dic, fit_inla_spatial_gs_cell$cpu.used[4],
fit_inla_spatial_vs_cell$waic$waic, fit_inla_spatial_vs_cell$dic$dic, fit_inla_spatial_vs_cell$cpu.used[4],
NA, NA, time_msgwr),
ncol = 3, byrow = TRUE))
names(df_comp) <- c("WAIC", "DIC", "Running Time")
row.names(df_comp) <- c("eq/stat terms",
"eq terms",
"no random effects",
"spatial eq/stat constants",
"spatial eq constant",
"spatial stat constants",
"spatial eq/stat constants + gs",
"spatial eq/stat constants + r (midpoint)",
"spatial eq/stat constants + r (eq coord)",
"spatial eq/stat constants + vs",
"full spatial model",
"spatial model, no random effects/constants",
"spatial eq/stat constants + cells",
"spatial eq/stat constants + cells, integration",
"spatial eq/stat constants + gs + vs + cells",
"spatial eq/stat constants + gs + cells",
"spatial eq/stat constants + vs + cells",
"MS-GWR")
knitr::kable(df_comp, digits = 2, row.names = TRUE,
caption = "Information Criteria (WAIC and DIC) for different INLA models
estimated on the Italian data.")| WAIC | DIC | Running Time | |
|---|---|---|---|
| eq/stat terms | -777.41 | -798.16 | 53.17 |
| eq terms | 2398.99 | 2392.42 | 42.54 |
| no random effects | 3407.34 | 3407.69 | 10.96 |
| spatial eq/stat constants | -1527.99 | -1511.66 | 805.22 |
| spatial eq constant | -1502.30 | -1473.62 | 157.87 |
| spatial stat constants | -1514.00 | -1497.32 | 184.21 |
| spatial eq/stat constants + gs | -1759.17 | -1736.45 | 800.53 |
| spatial eq/stat constants + r (midpoint) | -1627.22 | -1594.29 | 472.81 |
| spatial eq/stat constants + r (eq coord) | -1382.87 | -1427.63 | 401.44 |
| spatial eq/stat constants + vs | -1530.64 | -1512.75 | 1531.59 |
| full spatial model | -1780.31 | -1742.16 | 2131.77 |
| spatial model, no random effects/constants | -187.23 | -83.95 | 505.28 |
| spatial eq/stat constants + cells | -1841.46 | -1788.00 | 489.38 |
| spatial eq/stat constants + cells, integration | -1835.83 | -1784.71 | 2048.90 |
| spatial eq/stat constants + gs + vs + cells | -1964.89 | -1914.46 | 2152.52 |
| spatial eq/stat constants + gs + cells | -1955.72 | -1904.12 | 378.04 |
| spatial eq/stat constants + vs + cells | -1849.08 | -1801.96 | 951.04 |
| MS-GWR | NA | NA | 14189.79 |
All spatial models are an improvement over the simple ergodic model that only includes event and station terms. The spatial model that only includes spatially varying geometrical spreading, \(V_{S30}\)-scaling, ad linear R-scaling performs than the ergodic random effects model. The model using the spatially varying linear R term based on midpoints between events/stations performs better than the one using event coordinates; however, the models with cell-specific linear R coefficients outperform the spatially varying models.
8.2 Posterior distributions
posterior_inla <- as.data.frame(fit_inla$marginals.fixed$R)
posterior_inla_noranef <- as.data.frame(fit_inla_noranef$marginals.fixed$R)
posterior_inla_full <- as.data.frame(fit_inla_spatial_full$marginals.fixed$R)
posterior_inla_full_b <- as.data.frame(fit_inla_spatial_full_b$marginals.fixed$R)
posterior_inla_cell <- as.data.frame(fit_inla_spatial_cell$marginals.fixed$R)
posterior_inla_cell_int <- as.data.frame(fit_inla_spatial_cell_int$marginals.fixed$R)
posterior_inla_full_cell <- as.data.frame(fit_inla_spatial_full_cell$marginals.fixed$R)
ggplot() +
geom_line(posterior_inla, mapping = aes(x,y)) +
geom_line(posterior_inla_noranef, mapping = aes(x,y), colour = "red") +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_full_b, mapping = aes(x,y), colour = "cornflowerblue") +
geom_line(posterior_inla_cell, mapping = aes(x,y), colour = "cyan") +
geom_line(posterior_inla_cell_int, mapping = aes(x,y), colour = "cyan", linetype = 'dashed') +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "phi_SS", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of linear R coefficient.
posterior_inla <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla$internal.marginals.hyperpar[[1]]))
posterior_inla_noranef <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_noranef$internal.marginals.hyperpar[[1]]))
posterior_inla_full <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full$internal.marginals.hyperpar[[1]]))
posterior_inla_full_b <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full_b$internal.marginals.hyperpar[[1]]))
posterior_inla_cell <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell$internal.marginals.hyperpar[[1]]))
posterior_inla_cell_int <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell_int$internal.marginals.hyperpar[[1]]))
posterior_inla_full_cell <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full_cell$internal.marginals.hyperpar[[1]]))
ggplot() +
geom_line(posterior_inla, mapping = aes(x,y)) +
geom_line(posterior_inla_noranef, mapping = aes(x,y), colour = "red") +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_full_b, mapping = aes(x,y), colour = "cornflowerblue") +
geom_line(posterior_inla_cell, mapping = aes(x,y), colour = "cyan") +
geom_line(posterior_inla_cell_int, mapping = aes(x,y), colour = "cyan", linetype = 'dashed') +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "phi_SS", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of phi for different models.
posterior_inla <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla$internal.marginals.hyperpar[[2]]))
posterior_inla_full <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full$internal.marginals.hyperpar[[2]]))
posterior_inla_cell <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell$internal.marginals.hyperpar[[2]]))
posterior_inla_cell_int <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell_int$internal.marginals.hyperpar[[2]]))
posterior_inla_full_cell <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full_cell$internal.marginals.hyperpar[[2]]))
ggplot() +
geom_line(posterior_inla, mapping = aes(x,y)) +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_cell, mapping = aes(x,y), colour = "cyan") +
geom_line(posterior_inla_cell_int, mapping = aes(x,y), colour = "cyan", linetype = 'dashed') +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "tau", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of tau for different models.
posterior_inla <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla$internal.marginals.hyperpar[[3]]))
posterior_inla_full <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full$internal.marginals.hyperpar[[3]]))
posterior_inla_cell <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell$internal.marginals.hyperpar[[3]]))
posterior_inla_cell_int <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell_int$internal.marginals.hyperpar[[3]]))
posterior_inla_full_cell <- as.data.frame(inla.tmarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full_cell$internal.marginals.hyperpar[[3]]))
ggplot() +
geom_line(posterior_inla, mapping = aes(x,y)) +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_cell, mapping = aes(x,y), colour = "cyan") +
geom_line(posterior_inla_cell_int, mapping = aes(x,y), colour = "cyan", linetype = 'dashed') +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "phi_S2S", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of phi_S2S for different models.
posterior_inla_full <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full$internal.marginals.hyperpar$`log(Range) for idx.eq`))
posterior_inla_cell <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_cell$internal.marginals.hyperpar$`log(Range) for idx.eq`))
posterior_inla_cell_int <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_cell_int$internal.marginals.hyperpar$`log(Range) for idx.eq`))
posterior_inla_full_cell <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full_cell$internal.marginals.hyperpar$`log(Range) for idx.eq`))
ggplot() +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_cell, mapping = aes(x,y), colour = "cyan") +
geom_line(posterior_inla_cell_int, mapping = aes(x,y), colour = "cyan", linetype = 'dashed') +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "spatial range for event terms", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of spatial range of event constant for different models.
posterior_inla_full <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full$internal.marginals.hyperpar$`log(Range) for idx.stat`))
posterior_inla_cell <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_cell$internal.marginals.hyperpar$`log(Range) for idx.stat`))
posterior_inla_cell_int <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_cell_int$internal.marginals.hyperpar$`log(Range) for idx.stat`))
posterior_inla_full_cell <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full_cell$internal.marginals.hyperpar$`log(Range) for idx.stat`))
ggplot() +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_cell, mapping = aes(x,y), colour = "cyan") +
geom_line(posterior_inla_cell_int, mapping = aes(x,y), colour = "cyan", linetype = 'dashed') +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "spatial range for station terms", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of spatial range of station constant for different models.
posterior_inla_full <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full$internal.marginals.hyperpar$`log(Range) for idx.lnr`))
posterior_inla_full_b <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full_b$internal.marginals.hyperpar$`log(Range) for idx.lnr`))
posterior_inla_full_cell <- as.data.frame(inla.tmarginal(function(x) exp(x),
fit_inla_spatial_full_cell$internal.marginals.hyperpar$`log(Range) for idx.lnr`))
ggplot() +
geom_line(posterior_inla_full, mapping = aes(x,y), colour = "blue") +
geom_line(posterior_inla_full_b, mapping = aes(x,y), colour = "cornflowerblue") +
geom_line(posterior_inla_full_cell, mapping = aes(x,y), colour = "magenta") +
labs(x = "spatial range for geometrical spreading terms", y = "posterior") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Posterior distribution of spatial range of geometrical spreading terms for different models.
8.3 Compason of Fixed Effects
###########
# fixed effects plots
res <- fit_inla$summary.fixed[,c(4,3,5)]
res <- rbind(res, fit_inla_eq$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial_full$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial_full_cell$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial_cell$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial_full_b$summary.fixed[,c(4,3,5)])
colnames(res) = c("E", "L", "U")
rownames(res)=NULL
n.covar = nrow(fit_inla$summary.fixed)
model_list <- c("M1", "M2", "M3", "M4", "M5", "M6", "M7")
res$model = factor(rep(model_list, each=n.covar),
levels = model_list)
covar_list <- rownames(fit_inla$summary.fixed)
covar_list <- c("a", "b[1]", "b[2]", "c[2]", "c[1]", "c[3]", "f[1]", "f[2]", "k")
res$covar = factor(rep(covar_list, length(model_list)),
levels = covar_list)
ggplot(res, aes(x = model, y = E)) +
facet_wrap(~covar, scales = "free_y", labeller = label_parsed) +
geom_point(size = 3) +
geom_pointrange(aes(ymax = U, ymin = L)) +
xlab(NULL) + ylab(NULL) +
theme(
axis.title = element_text(size = 18),
axis.text = element_text(size = 12),
plot.title = element_text(size = 30),
strip.text.x = element_text(size = 12)
)
Comparison of fixed effects.
8.4 Comparison of Standard Deviations
# calculate expected value of standard deviations from posterior marginals of log precisions
tmp <- matrix(c(inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla$internal.marginals.hyperpar$`Log precision for eq`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla$internal.marginals.hyperpar$`Log precision for stat`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_noranef$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
NA, NA,
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial$internal.marginals.hyperpar$`Log precision for eq`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial$internal.marginals.hyperpar$`Log precision for stat`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full$internal.marginals.hyperpar$`Log precision for eq`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full$internal.marginals.hyperpar$`Log precision for stat`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_full_b$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
NA, NA,
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell$internal.marginals.hyperpar$`Log precision for eq`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell$internal.marginals.hyperpar$`Log precision for stat`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell_int$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell_int$internal.marginals.hyperpar$`Log precision for eq`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_spatial_cell_int$internal.marginals.hyperpar$`Log precision for stat`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_resid$internal.marginals.hyperpar$`Log precision for the Gaussian observations`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_resid$internal.marginals.hyperpar$`Log precision for eq`),
inla.emarginal(function(x) sqrt(exp(-x)),
fit_inla_resid$internal.marginals.hyperpar$`Log precision for stat`),
0.220582, 0.155988, 0.200099
), ncol = 3, byrow = TRUE)
row.names(tmp) <- c("eq/stat terms",
"no random effects",
"spatial eq/stat constants",
"full spatial model",
"spatial model, no random effects/constants",
"spatial eq/stat constants + cells",
"spatial eq/stat constants + cells, integration",
"MS-GWR residuals",
"ITA18"
)
knitr::kable(tmp, col.names = c("phi_SS","tau","phi_S2S"), row.names = TRUE,
caption = "Estimated standard deviations (mean of posterior distribution).")| phi_SS | tau | phi_S2S | |
|---|---|---|---|
| eq/stat terms | 0.2042613 | 0.1449961 | 0.2331139 |
| no random effects | 0.3450781 | NA | NA |
| spatial eq/stat constants | 0.1821196 | 0.0829741 | 0.1590015 |
| full spatial model | 0.1770333 | 0.0831417 | 0.1545714 |
| spatial model, no random effects/constants | 0.2194835 | NA | NA |
| spatial eq/stat constants + cells | 0.1751112 | 0.0862417 | 0.1573863 |
| spatial eq/stat constants + cells, integration | 0.1749622 | 0.0873041 | 0.1570731 |
| MS-GWR residuals | 0.2003396 | 0.0916589 | 0.2066721 |
| ITA18 | 0.2205820 | 0.1559880 | 0.2000990 |
8.5 Model Comparison based on Cross-Validation
Now we calculate a 5-fold CV for some selected models (5 folds because of time). The CV for the MS-GWR is calculated outside, since it takes a long time (but is run with the same folds). We randomize the event indices, and partition them into five sets. We loop over sets, and for each set, set the value of the response for the test set to NA. inla will automatically calculate predictions for data where the response is NA. We then calculate the residuals for the test set, and store the test residuals. We also calculate RMSE for each fold.
For GMMs, is important that the variability is well estimated, no just median predictions. Hence, it is important to not only compare es residuals, bu also the likelihood of the test data given the model. Here, we calculate the loglikelihood of the test data, and sample from the posterior distributions of the hyperparameters, and average over the samples for each test record. This part is only done for the INLA models.
library(matrixStats) # for logsumexp
# load cv results of MS-GWR
load(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY', 'RESULTS', 'results_cveq_fold_b_ll_mswgr.Rdata'))
tmp <- resid_mod
tmp2 <- cv_results
tmp3 <- loglik_mod
n_eq <- max(data_reg$eq)
set.seed(1701)
idx_rand <- sample(1:n_eq)
n_fold <- 5
batch = floor(n_eq/n_fold)
rest = n_eq - batch * n_fold
rest_const = rest
n_sample <- 1000
n_mod <- 10
cv_results <- matrix(nrow = n_mod, ncol = n_fold)
lpd_results <- matrix(nrow = n_mod, ncol = n_fold)
resid_mod <- matrix(nrow = n_rec, ncol = n_mod)
loglik_mod <- matrix(nrow = n_rec, ncol = n_mod)
# last two models are MS-GWR
resid_mod[,c(n_mod - 1, n_mod)] <- tmp
loglik_mod[,c(n_mod - 1, n_mod)] <- tmp3
cv_results[c(n_mod - 1, n_mod),] <- tmp2
for(b in 1:n_fold) {
print(paste0('fold ', b))
if (rest >0){
batch_t = batch + 1
beginning = batch_t*(b-1)+1
end = b*batch_t
rest = rest - 1
} else {
beginning = rest_const + batch*(b-1) + 1
end = rest_const + batch*b
}
idx_test_eq <- idx_rand[beginning:end]
idx_test <- which(data_reg$eq %in% idx_test_eq)
y <- data_reg$Y
y[idx_test] <- NA
data_reg_p <- data_reg
data_reg_p$Y <- y
### ergodic model with random effects
mod_idx <- 1
print(paste0('working on model ', mod_idx))
form_p <- Y ~ M1 + M2 + lnR + MlnR + R + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg)
fit_inla_p <- inla(form_p,
data = data_reg_p,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
control.predictor = list(compute = TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### ergodic model with random effects for events
mod_idx <- 2
print(paste0('working on model ', mod_idx))
form_p <- Y ~ M1 + M2 + lnR + MlnR + R + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg)
fit_inla_p <- inla(form_p,
data = data_reg_p,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
control.predictor = list(compute = TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### ergodic model without random effects
mod_idx <- 3
print(paste0('working on model ', mod_idx))
fit_inla_p <- inla(Y ~ M1 + M2 + lnR + MlnR + R + lnVS,
data = data_reg_p,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
control.predictor = list(compute = TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### spatial model without random effects
mod_idx <- 4
print(paste0('working on model ', mod_idx))
theta <- fit_inla_spatial_full_b$mode$theta
stk_spatial_p <- inla.stack(
data = list(y = y),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_lnr, A_r, A_vs, 1),
effects = list(idx.lnr = idx.lnr,
idx.r = idx.r,
idx.vs = idx.vs,
data_reg_p
),
tag = 'model_spatial_p')
form_p <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + lnVS +
f(idx.lnr, model = spde_lnR) + f(idx.r, model = spde_R) +
f(idx.vs, model = spde_vs)
fit_inla_p <- inla(form_p,
data = inla.stack.data(stk_spatial_p),
control.predictor = list(A = inla.stack.A(stk_spatial_p),
compute = TRUE),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### spatial model with event/station constants
mod_idx <- 5
print(paste0('working on model ', mod_idx))
theta <- fit_inla_spatial$mode$theta
stk_spatial_p <- inla.stack(
data = list(y = y),
A = list(A_eq, A_stat, 1),
effects = list(idx.eq = idx.eq,
idx.stat = idx.stat,
data_reg_p
),
tag = 'model_spatial_p')
form_p <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq)
fit_inla_p <- inla(form_p,
data = inla.stack.data(stk_spatial_p),
control.predictor = list(A = inla.stack.A(stk_spatial_p),
compute = TRUE),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### full spatial model with cells
mod_idx <- 6
print(paste0('working on model ', mod_idx))
theta <- fit_inla_spatial_full_cell$mode$theta
stk_spatial_p <- inla.stack(
data = list(y = y),
A = list(A_eq, A_stat, A_lnr, A_vs, 1),
effects = list(idx.eq = idx.eq,
idx.stat = idx.stat,
idx.lnr = idx.lnr,
idx.vs = idx.vs,
data_reg_p
),
tag = 'model_spatial_p')
form_p <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.lnr, model = spde_lnR) + f(idx.vs, model = spde_vs) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
fit_inla_p <- inla(form_p,
data = inla.stack.data(stk_spatial_p),
control.predictor = list(A = inla.stack.A(stk_spatial_p),
compute = TRUE),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### spatial model with eq/stat, gs and cells
mod_idx <- 7
print(paste0('working on model ', mod_idx))
theta <- fit_inla_spatial_gs_cell$mode$theta
stk_spatial_p <- inla.stack(
data = list(y = y),
A = list(A_eq, A_stat, A_lnr, 1),
effects = list(idx.eq = idx.eq,
idx.stat = idx.stat,
idx.lnr = idx.lnr,
data_reg_p
),
tag = 'model_spatial_p')
form_p <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx.lnr, model = spde_lnR) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
fit_inla_p <- inla(form_p,
data = inla.stack.data(stk_spatial_p),
control.predictor = list(A = inla.stack.A(stk_spatial_p),
compute = TRUE),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
### spatial model with constants and cells
mod_idx <- 8
print(paste0('working on model ', mod_idx))
theta <- fit_inla_spatial_cell$mode$theta
stk_spatial_p <- inla.stack(
data = list(y = y),
A = list(A_eq, A_stat, 1),
effects = list(idx.eq = idx.eq,
idx.stat = idx.stat,
data_reg_p
),
tag = 'model_spatial_p')
form_p <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx_cell, model = "z", Z = dm_sparse, hyper = prior_prec_cell)
fit_inla_p <- inla(form_p,
data = inla.stack.data(stk_spatial_p),
control.predictor = list(A = inla.stack.A(stk_spatial_p),
compute = TRUE),
control.compute = list(dic = TRUE, cpo = TRUE, waic = TRUE, config = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian"),
control.mode=list(theta = theta, restart=TRUE)
)
pred <- fit_inla_p$summary.fitted.values[idx_test,]$mean
y_res <- data_reg$Y[idx_test] - pred
cv_results[mod_idx, b] <- sqrt(mean(y_res^2))
resid_mod[idx_test, mod_idx] <- y_res
## sample
sample <- inla.posterior.sample(n_sample, fit_inla_p)
y_sim <- matrix(ncol = length(idx_test), nrow = n_sample)
tau_sim <- vector(length = n_sample)
for(i in 1:n_sample) {
y_sim[i,] <- sample[[i]]$latent[idx_test]
tau_sim[i] <- sample[[i]]$hyperpar[1]
}
residuals <- sweep(y_sim, 2, data_reg$Y[idx_test]) #subtract true values from each row of y_sim
log_dens <- dnorm(residuals, sd=sqrt(1/tau_sim), log=TRUE) #log likelihood of the residual of each draw
lpd <- apply(log_dens, 2, function (col) {logSumExp(col) - log(length(col))}) #take mean across draws (rows) for each point, giving pointwise likelihoods; using log_sum_exp in case of underflow
loglik_mod[idx_test, mod_idx] <- lpd
lpd_results[mod_idx, b] <- sum(lpd)
print(paste0('mean = ',mean(y_res), ' sd = ',sd(y_res)))
rm(fit_inla_p)
}## [1] "fold 1"
## [1] "working on model 1"
## [1] "mean = 0.00456364312555309 sd = 0.246874031197999"
## [1] "working on model 2"
## [1] "mean = -0.020944816580942 sd = 0.322897634952355"
## [1] "working on model 3"
## [1] "mean = 0.0120734080940625 sd = 0.325950242557458"
## [1] "working on model 4"
## [1] "mean = -0.00576785088628497 sd = 0.290883912877798"
## [1] "working on model 5"
## [1] "mean = 0.0105704062171967 sd = 0.240168217470607"
## [1] "working on model 6"
## [1] "mean = 0.0270167645024283 sd = 0.244133208817596"
## [1] "working on model 7"
## [1] "mean = 0.0254735889877114 sd = 0.24568415697194"
## [1] "working on model 8"
## [1] "mean = 0.00955557725266701 sd = 0.239121190719968"
## [1] "fold 2"
## [1] "working on model 1"
## [1] "mean = 0.00404755503513967 sd = 0.280325067852517"
## [1] "working on model 2"
## [1] "mean = -0.0117057507940665 sd = 0.358996694938169"
## [1] "working on model 3"
## [1] "mean = 0.00209357611810311 sd = 0.359328553663907"
## [1] "working on model 4"
## [1] "mean = 0.0444849274594708 sd = 0.278249026951578"
## [1] "working on model 5"
## [1] "mean = 0.0155597344999331 sd = 0.27111527631352"
## [1] "working on model 6"
## [1] "mean = 0.0373896509679996 sd = 0.26470761088567"
## [1] "working on model 7"
## [1] "mean = 0.0337482959798949 sd = 0.263891483852901"
## [1] "working on model 8"
## [1] "mean = 0.0125158229191847 sd = 0.265052788042987"
## [1] "fold 3"
## [1] "working on model 1"
## [1] "mean = 0.0381252043879595 sd = 0.28736787487984"
## [1] "working on model 2"
## [1] "mean = 0.0599004736792779 sd = 0.364538069524996"
## [1] "working on model 3"
## [1] "mean = 0.0898962404088453 sd = 0.367356927845064"
## [1] "working on model 4"
## [1] "mean = -0.04829865021598 sd = 0.340966218331341"
## [1] "working on model 5"
## [1] "mean = 0.0385085278731682 sd = 0.270522254554355"
## [1] "working on model 6"
## [1] "mean = 0.0362427153629523 sd = 0.260967891746895"
## [1] "working on model 7"
## [1] "mean = 0.0303657781372568 sd = 0.262978485034134"
## [1] "working on model 8"
## [1] "mean = 0.0412424605700754 sd = 0.262546668885539"
## [1] "fold 4"
## [1] "working on model 1"
## [1] "mean = -0.0413461516430655 sd = 0.245148863948091"
## [1] "working on model 2"
## [1] "mean = -0.0795669678083998 sd = 0.322666503334188"
## [1] "working on model 3"
## [1] "mean = -0.0640222761949934 sd = 0.324012333043252"
## [1] "working on model 4"
## [1] "mean = 0.0402718561206015 sd = 0.321990631952038"
## [1] "working on model 5"
## [1] "mean = -0.00764105879713926 sd = 0.245321295376291"
## [1] "working on model 6"
## [1] "mean = 0.00854846079373443 sd = 0.26861900142539"
## [1] "working on model 7"
## [1] "mean = 0.00859157328779772 sd = 0.268194869748081"
## [1] "working on model 8"
## [1] "mean = 0.00548000013469638 sd = 0.244772094682707"
## [1] "fold 5"
## [1] "working on model 1"
## [1] "mean = -0.0775772475564917 sd = 0.28792767502906"
## [1] "working on model 2"
## [1] "mean = -0.0874268291506236 sd = 0.366153595890512"
## [1] "working on model 3"
## [1] "mean = -0.0759052356549884 sd = 0.366376078721184"
## [1] "working on model 4"
## [1] "mean = -0.0261302059107283 sd = 0.268980920484973"
## [1] "working on model 5"
## [1] "mean = -0.0462927673222794 sd = 0.253871830836452"
## [1] "working on model 6"
## [1] "mean = -0.0433563120966332 sd = 0.24299059021757"
## [1] "working on model 7"
## [1] "mean = -0.0426394338411577 sd = 0.243832712348802"
## [1] "working on model 8"
## [1] "mean = -0.0475885460833465 sd = 0.251452781712103"row.names(cv_results) <- c("eq/stat random effects", "eq random effects", "no random effects",
"spatial without random effects", "spatial eq/stat constants",
"full spatial model + cells", "eq/stat constants + gs + cells",
"eq/stat constants + cells", "MS-GWR", "MS-GWR stat")
row.names(lpd_results) <- row.names(cv_results)
knitr::kable(cv_results, col.names = c("fold 1", "fold 2", "fold 3", "fold 4", "fold 5"),
caption = "Root mean square error (RMSE) for individual folds from cross-validation.")| fold 1 | fold 2 | fold 3 | fold 4 | fold 5 | |
|---|---|---|---|---|---|
| eq/stat random effects | 0.2468171 | 0.2801823 | 0.2897683 | 0.2484781 | 0.2979645 |
| eq random effects | 0.3234469 | 0.3589673 | 0.3692783 | 0.3321597 | 0.3761505 |
| no random effects | 0.3260430 | 0.3591141 | 0.3780491 | 0.3301020 | 0.3738583 |
| spatial without random effects | 0.2908244 | 0.2816140 | 0.3442307 | 0.3243235 | 0.2700247 |
| spatial eq/stat constants | 0.2403044 | 0.2713953 | 0.2731388 | 0.2453054 | 0.2578505 |
| full spatial model + cells | 0.2455262 | 0.2671743 | 0.2633659 | 0.2686073 | 0.2466295 |
| eq/stat constants + gs + cells | 0.2469031 | 0.2658801 | 0.2646180 | 0.2681850 | 0.2473333 |
| eq/stat constants + cells | 0.2392161 | 0.2651856 | 0.2656592 | 0.2446988 | 0.2557110 |
| MS-GWR | 0.3027971 | 0.3388769 | 0.3423976 | 0.3019265 | 0.3428962 |
| MS-GWR stat | 0.2530393 | 0.2695278 | 0.2763125 | 0.2449033 | 0.2972949 |
knitr::kable(lpd_results, col.names = c("fold 1", "fold 2", "fold 3", "fold 4", "fold 5"),
caption = "Log-likelihood for individual folds from cross-validation.")| fold 1 | fold 2 | fold 3 | fold 4 | fold 5 | |
|---|---|---|---|---|---|
| eq/stat random effects | -40.00602 | -110.74418 | -220.4781 | -26.509153 | -122.775796 |
| eq random effects | -377.10096 | -324.55453 | -516.2534 | -296.297332 | -267.638763 |
| no random effects | -379.79510 | -323.37943 | -557.8988 | -285.619707 | -267.245925 |
| spatial without random effects | -203.20551 | -167.62366 | -477.3080 | -297.398313 | -84.165858 |
| spatial eq/stat constants | 27.99359 | -81.27059 | -149.5843 | 6.708486 | -38.380161 |
| full spatial model + cells | 10.60474 | -94.39371 | -197.5892 | -61.539127 | -6.687196 |
| eq/stat constants + gs + cells | 4.14287 | -87.32077 | -202.3223 | -58.540895 | -7.669160 |
| eq/stat constants + cells | 42.36296 | -60.58785 | -142.5026 | 17.413093 | -34.209705 |
| MS-GWR | NA | NA | NA | NA | NA |
| MS-GWR stat | NA | NA | NA | NA | NA |
res <- t(rbind(sqrt(colMeans(resid_mod[,]^2)),
colSums(loglik_mod)))
row.names(res) <- row.names(cv_results)
knitr::kable(res, col.names = c("RMSE", "LL"),
caption = "RMSE and total loglikelihood on test set from cross-validation.")| RMSE | LL | |
|---|---|---|
| eq/stat random effects | 0.2709532 | -520.5132 |
| eq random effects | 0.3500225 | -1781.8450 |
| no random effects | 0.3523609 | -1813.9390 |
| spatial without random effects | 0.3077671 | -1229.7014 |
| spatial eq/stat constants | 0.2574662 | -234.5330 |
| full spatial model + cells | 0.2584512 | -349.6045 |
| eq/stat constants + gs + cells | 0.2588904 | -351.7103 |
| eq/stat constants + cells | 0.2537116 | -177.5241 |
| MS-GWR | 0.3244398 | -1420.5269 |
| MS-GWR stat | 0.2662948 | NA |
9 Comparisons of Spatial Coefficients
### grid from mwgr
co_df_utm <- coords_df_utm/1000
proj <- inla.mesh.projector(mesh, loc = as.matrix(co_df_utm))9.1 Plot of Spatially Varying Event Terms
field.proj <- inla.mesh.project(proj, fit_inla_spatial_cell$summary.random$idx.eq$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_c <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.2,0.2), breaks = c(-0.2, -0.1,0, 0.1, 0.2)
) +
labs(title = "Spatial event term - mean",
subtitle = "spatial eq/stat + cell model",
x = "UTM X (km)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_cell$summary.random$idx.eq$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_c <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.25), breaks = c(0, 0.1, 0.2)
) +
labs(title = "Spatial event term - sd",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.eq$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.2,0.2), breaks = c(-0.2, -0.1,0, 0.1, 0.2)
) +
labs(title = "Spatial event term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.eq$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.25), breaks = c(0, 0.1, 0.2)
) +
labs(title = "Spatial event term - sd",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_c, p_sd_c, p_mean_f, p_sd_f, nrow = 2)
9.2 Plot of Spatially Varying Station Terms
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.stat$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.5,0.5), breaks = c(-0.5, -0.25,0, 0.25, 0.5)
) +
labs(title = "Spatial station term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.stat$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.25), breaks = c(0, 0.1, 0.2)
) +
labs(title = "Spatial station term - sd",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
9.3 Plot of Spatially Varying Geometrical Spreading
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.lnr$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.2,0.2), breaks = c(-0.2, -0.1,0, 0.1, 0.2)
) +
labs(title = "Spatial GS term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.lnr$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.2), breaks = c(0, 0.1, 0.2)
) +
labs(title = "Spatial GS term - sd",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# spatial without random effects and constants
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.lnr$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.2,0.2), breaks = c(-0.2, -0.1,0, 0.1, 0.2)
) +
labs(title = "Spatial delta GS term - mean",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.lnr$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.2), breaks = c(0, 0.1, 0.2)
) +
labs(title = "Spatial delta GS term - sd",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_f, p_sd_f, p_mean_fb, p_sd_fb, nrow = 2)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.lnr$mean +
fit_inla_spatial_full$summary.fixed$mean[4])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(title = "Spatial GS term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# spatial model without random effects and constants
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.lnr$mean +
fit_inla_spatial_full_b$summary.fixed$mean[4])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(title = "Spatial GS term - mean",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# full spatial model with cells
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_cell$summary.random$idx.lnr$mean +
fit_inla_spatial_full_cell$summary.fixed$mean[4])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fc <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(title = "Spatial GS term - mean",
subtitle = "full spatial model with cells",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# MS-WGR model
df_plot <- as.data.frame(cbind(co_df_utm, beta_c2))
p_mean_wgr <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=beta_c2))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(title = "Spatial GS term",
subtitle = "MS-WGR model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_f, p_mean_fb, p_mean_fc, p_mean_wgr, nrow = 2)
9.4 Plot of Spatially Varying Anelastic Attenuation
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.r$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.004,0.004), breaks = c(-0.004, -0.002,0, 0.002, 0.004)
) +
labs(title = "Spatial delta linear R term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.r$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.003), breaks = c(0, 0.0015, 0.003)
) +
labs(title = "Spatial delta linear R term - sd",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# spatial without random effects and constants
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.r$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.004,0.004), breaks = c(-0.004, -0.002,0, 0.002, 0.004)
) +
labs(title = "Spatial delta linear R term - mean",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.r$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.003), breaks = c(0, 0.0015, 0.003)
) +
labs(title = "Spatial delta linear R term - sd",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# cell-specific attenuation
cells <- read.csv(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/',
'DATA', 'cells_dim25.csv'))
cell_idx <- data_dm$cell_idx
cells2 <- cells[cell_idx,]
ca <- fit_inla_spatial_full_cell$summary.random$idx_cell$mean[(n_rec + 1):(n_rec + n_cell)]
ca_sd <- fit_inla_spatial_full_cell$summary.random$idx_cell$sd[(n_rec + 1):(n_rec + n_cell)]
df_plot <- data.frame(x = cells2$X0, y = cells2$X1,
z = ca, z_sd = ca_sd)
p_mean_c <- ggplot() +
geom_raster(df_plot, mapping = aes(x=x, y=y, fill=z))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "cells"
,limits = c(-0.004,0.004), breaks = c(-0.004, -0.002,0, 0.002, 0.004)
) +
labs(title = "Cell-specific linear R term - mean",
subtitle = "full spatial model with cells",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
p_sd_c <- ggplot() +
geom_raster(df_plot, mapping = aes(x=x, y=y, fill=z_sd))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "cells"
,limits = c(0,0.003), breaks = c(0, 0.0015, 0.003)
) +
labs(title = "Cell-specific linear R term - sd",
subtitle = "full spatial model with cells",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_f, p_sd_f, p_mean_fb, p_sd_fb, p_mean_c, p_sd_c, nrow = 3)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.r$mean +
fit_inla_spatial_full$summary.fixed$mean[6])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.009,0.003), breaks = c(-0.009, -0.004, 0.003)
) +
labs(title = "Spatial linear term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# spatial model without random effects and constants
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.r$mean +
fit_inla_spatial_full_b$summary.fixed$mean[6])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.009,0.003), breaks = c(-0.009, -0.004, 0.003)
) +
labs(title = "Spatial linear term - mean",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# MS-WGR model
df_plot <- as.data.frame(cbind(co_df_utm, beta_c3))
p_mean_wgr <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=beta_c3))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.009,0.003), breaks = c(-0.009, -0.004, 0.003)
) +
labs(title = "Spatial linear term",
subtitle = "MS-WGR model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# cell-specific model
df_plot <- data.frame(x = cells2$X0, y = cells2$X1,
z = ca + fit_inla_spatial_full_cell$summary.fixed$mean[6])
p_mean_c <- ggplot() +
geom_raster(df_plot, mapping = aes(x=x, y=y, fill=z))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "cells"
,limits = c(-0.009,0.003), breaks = c(-0.009, -0.004, 0.003)
) +
labs(title = "Cell-specific linear R term ",
subtitle = "full spatial model with cells",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_f, p_mean_fb, p_mean_wgr, p_mean_c, nrow = 2)
9.5 Plot of Spatially Varying VS30 Scaling Term
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.vs$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.6,0.6), breaks = c(-0.6, -0.3,0, 0.3, 0.6)
) +
labs(title = "Spatial delta VS - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.vs$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.6), breaks = c(0, 0.3, 0.6)
) +
labs(title = "Spatial delta VS term - sd",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# spatial without random effects and constants
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.vs$mean)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-0.6,0.6), breaks = c(-0.6, -0.3,0, 0.3, 0.6)
) +
labs(title = "Spatial delta VS term - mean",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.vs$sd)
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_sd_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(0,0.6), breaks = c(0, 0.3, 0.6)
) +
labs(title = "Spatial delta VS term - sd",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_f, p_sd_f, p_mean_fb, p_sd_fb, nrow = 2)
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.vs$mean +
fit_inla_spatial_full$summary.fixed$mean[9])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_f <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.,0.00), breaks = c(-1., -0.5, 0.)
) +
labs(title = "Spatial VS term - mean",
subtitle = "full spatial model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# spatial model without random effects and constants
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.vs$mean +
fit_inla_spatial_full_b$summary.fixed$mean[9])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p_mean_fb <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.,0.00), breaks = c(-1., -0.5, 0.)
) +
labs(title = "Spatial VS term - mean",
subtitle = "full spatial model without random effects",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
# MS-WGR model
df_plot <- as.data.frame(cbind(co_df_utm, beta_k))
p_mean_wgr <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=beta_k))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = ""
,limits = c(-1.,0.00), breaks = c(-1., -0.5, 0.)
) +
labs(title = "Spatial VS term",
subtitle = "MS-WGR model",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_f, p_mean_fb, p_mean_wgr, nrow = 2)
9.6 Comparison of Cell-Specific Attenuation models
ca <- fit_inla_spatial_full_cell$summary.random$idx_cell$mean[(n_rec + 1):(n_rec + n_cell)]
ca_sd <- fit_inla_spatial_full_cell$summary.random$idx_cell$sd[(n_rec + 1):(n_rec + n_cell)]
ca2 <- fit_inla_spatial_cell$summary.random$idx_cell$mean[(n_rec + 1):(n_rec + n_cell)]
ca2_sd <- fit_inla_spatial_cell$summary.random$idx_cell$sd[(n_rec + 1):(n_rec + n_cell)]
df_plot <- data.frame(x = cells2$X0, y = cells2$X1,
z = ca, z_sd = ca_sd,
z2 = ca2, z2_sd = ca2_sd,
z_diff = ca - ca2)
p_mean_c <- ggplot() +
geom_raster(df_plot, mapping = aes(x=x, y=y, fill=z))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "cells"
,limits = c(-0.004,0.004), breaks = c(-0.004, -0.002,0, 0.002, 0.004)
) +
labs(title = "Cell-specific linear R term - mean",
subtitle = "full spatial model with cells",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
p_mean_c2 <- ggplot() +
geom_raster(df_plot, mapping = aes(x=x, y=y, fill=z2))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "cells"
,limits = c(-0.004,0.004), breaks = c(-0.004, -0.002,0, 0.002, 0.004)
) +
labs(title = "Cell-specific linear R term - mean",
subtitle = "Spatial model with constants and cells",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
p_mean_diff <- ggplot() +
geom_raster(df_plot, mapping = aes(x=x, y=y, fill=z_diff))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "cells"
,limits = c(-0.004,0.004), breaks = c(-0.004, -0.002,0, 0.002, 0.004)
) +
labs(title = "Cell-specific linear R term - difference",
subtitle = "",
x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
grid.arrange(p_mean_c, p_mean_c2, p_mean_diff, nrow = 2)
10 Prediction
Here, we calculate predictions for some example scenarios, and compare different models. We can make predictions by refitting different models and adding the scenarios for which we want to predict, with the target variable \(y\) set to NA. In that case (and if we add control.predictor = list(compute = TRUE) to the INLA options), the posterior distribution of \(y\) is calculated.
Below, we load some data for which want to make predictions. There are three different events (North, Middle, South), with stations from 1km to 200km in different directions for each (East, North, West, South). The scenarios \(M_W = 6\), \(V_{S30} = 400\), \(F_{SS} = 1\). We set he event index to NA for the new data, which means that we do not include this term in the prediction. This means that \(\tau\) is no included in the predictive uncertainty.
data_pred <- read.csv("/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/data_prediction.csv")
data_dm_comb <- rstan::read_rdump('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/dm_combined_25x25.Rdata')
dm_sparse_c <- as(data_dm_comb$RC,"dgCMatrix")
load(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY','RESULTS', 'results_pred_mswgr.Rdata'))
attach(data_pred)
b1p = (mag-mh)*(mag<=mh)
b2p = (mag-mh)*(mag>mh)
c1p = (mag-mref)*log10(sqrt(JB_complete^2+h^2))
c2p = log10(sqrt(JB_complete^2+h^2))
c3p = sqrt(JB_complete^2+h^2)
f1p = as.numeric(fm_type_code == "SS")
f2p = as.numeric(fm_type_code == "TF")
kp = log10(vs30/800)*(vs30<=1500)+log10(1500/800)*(vs30>1500)
yp = log10(rotD50_pga)
detach(data_pred)
data_pred_reg <- data.frame(Y = yp,
M1 = b1p,
M2 = b2p,
MlnR = c1p,
lnR = c2p,
R = c3p,
Fss = f1p,
Frv = f2p,
lnVS = kp,
#eq = data_pred$EQID,
eq = NA,
stat = data_pred$STATID,
intercept = 1,
resid = NA
)
data_pred_reg$idx_cell <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
data_pred_reg$idx <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
co_eq_pred <- data_pred[,c(13,14)]
co_stat_pred <- data_pred[,c(15,16)]
co_eq_stat_pred <- cbind(rowMeans(cbind(co_eq_pred$X_ev, co_stat_pred$X_stat)),
rowMeans(cbind(co_eq_pred$Y_ev, co_stat_pred$Y_stat)))
ggplot() + theme_bw() + gg(mesh) +
geom_point(data = as.data.frame(co_stat_pred), aes(x=X_stat,y=Y_stat), color = set1[1]) +
geom_point(data = as.data.frame(co_eq_pred), aes(x=X_ev,y=Y_ev), color = set1[2]) +
labs(x="X (km)", y="Y (km)") +
theme(axis.title = element_text(size=30), axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20))
We alse calculate the \(A\) matrices which porject the spaial field from the mesh nodes to the prediction locations. Then we define the stack for the new locations, and combine the data and prediction stack.
A_stat_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_stat_pred))
A_eq_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_pred))
idx.stat_pred <- inla.spde.make.index("idx.stat",spde_stat$n.spde)
idx.eq_pred <- inla.spde.make.index("idx.eq",spde_eq$n.spde)
A_lnr_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_pred),
weights = data_pred_reg$lnR)
A_r_es_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_stat_pred),
weights = data_pred_reg$R)
idx.lnr_pred <- inla.spde.make.index("idx.lnr",spde_lnR$n.spde)
idx.r_pred <- inla.spde.make.index("idx.r",spde_R$n.spde)
A_vs_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_stat_pred),
weights = data_pred_reg$lnVS)
idx.vs_pred <- inla.spde.make.index("idx.vs",spde_vs$n.spde)
# create the stack
stk_spatial_full_pred <- inla.stack(
data = list(y = NA),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_eq_pred, A_lnr_pred, A_r_es_pred, A_stat_pred, A_vs_pred, 1),
effects = list(idx.eq = idx.eq_pred,
idx.lnr = idx.lnr_pred,
idx.r = idx.r_pred,
idx.stat = idx.stat_pred,
idx.vs = idx.vs_pred,
data_pred_reg
),
tag = 'model_spatial_full_pred')
stk_spatial_full_comb <- inla.stack(stk_spatial_full, stk_spatial_full_pred)Now we fit the ergodic model.
data_comb <- rbind(data_reg, data_pred_reg)
fit_inla_pred <- inla(form,
data = data_comb,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.predictor = list(compute = TRUE)
)
idx_pred <- inla.stack.index(stk_spatial_full_comb, 'model_spatial_full_pred')$data
pred <- fit_inla_pred$summary.fitted.values[idx_pred,]We fit the full spatial model, which includes all spatially varying coefficients (constants, geometrical spreading, linear R scaling, linear site scaling).
fit_inla_spatial_full_pred <- inla(form_spatial_full,
data = inla.stack.data(stk_spatial_full_comb),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)
idx_pred <- inla.stack.index(stk_spatial_full_comb, 'model_spatial_full_pred')$data
pred_spatial_full <- fit_inla_spatial_full_pred$summary.fitted.values[idx_pred,]Now we fit the spatial model with cell-specific attenuation. This model only includes event and station constants.
form_spatial_cell_c <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx_cell, model = "z", Z = dm_sparse_c, hyper = prior_prec_cell)
fit_inla_spatial_cell_pred <- inla(form_spatial_cell_c,
data = inla.stack.data(stk_spatial_full_comb),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'z' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.idx_pred <- inla.stack.index(stk_spatial_full_comb, 'model_spatial_full_pred')$data
pred_spatial_cell <- fit_inla_spatial_cell_pred$summary.fitted.values[idx_pred,]We fit the model again with predictions, but now we turn of he spatially varying station term, by setting the index to NA.
idx.stat_pred2 <- idx.stat_pred
idx.stat_pred2$idx.stat <- rep(NA, length(idx.stat_pred$idx.stat))
# create the stack
stk_spatial_full_pred2 <- inla.stack(
data = list(y = NA),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_eq_pred, A_lnr_pred, A_r_es_pred, A_stat_pred, A_vs_pred, 1),
effects = list(idx.eq = idx.eq_pred,
idx.lnr = idx.lnr_pred,
idx.r = idx.r_pred,
idx.stat = idx.stat_pred2,
idx.vs = idx.vs_pred,
data_pred_reg
),
tag = 'model_spatial_full_pred2')
stk_spatial_full_comb2 <- inla.stack(stk_spatial_full, stk_spatial_full_pred2)
fit_inla_spatial_cell_pred2 <- inla(form_spatial_cell_c,
data = inla.stack.data(stk_spatial_full_comb2),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb2),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'z' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.Now we plot the predictions. We plot the predictions for each event spearately, with different colors in for predictions in different directions.
df_plot <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, pred))
names(df_plot)[c(1,2,3,4)] <- c('EQID', 'WESN', 'R', 'Y_gwr')
df_plot_c <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, pred_spatial_cell))
names(df_plot_c)[c(1,2,3,4)] <- c('EQID', 'WESN', 'R', 'Y_gwr')
df_plot_f <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, pred_spatial_full))
names(df_plot_f)[c(1,2,3,4)] <- c('EQID', 'WESN', 'R', 'Y_gwr')
# # EQ 138
ref_id <- 138
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_f2 <- df_plot_f[df_plot_f$EQID == ref_id
& df_plot_f$WESN != 0,]
p138 <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
geom_line(df_plot_f2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN)), linetype = 'twodash') +
scale_x_continuous(trans='log10') +
ylim(-0.75, 3) +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
# EQ 139
ref_id <- 139
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_f2 <- df_plot_f[df_plot_f$EQID == ref_id
& df_plot_f$WESN != 0,]
p139 <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
geom_line(df_plot_f2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN)), linetype = 'twodash') +
scale_x_continuous(trans='log10') +
ylim(-0.75, 3) +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
# EQ 140
ref_id <- 140
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_f2 <- df_plot_f[df_plot_f$EQID == ref_id
& df_plot_f$WESN != 0,]
p140 <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
geom_line(df_plot_f2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN)), linetype = 'twodash') +
scale_x_continuous(trans='log10') +
ylim(-0.75, 3) +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
grid.arrange(p138, p139, p140, nrow = 2)
Now we plot the cell-specific model with and without the spatially varying station terms. This allows us to check which difference is due to the cell-specific attenuation, and what is due to station terms.
df_plot_c <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
pred_spatial_cell))
names(df_plot_c)[c(1,2,3)] <- c('EQID', 'WESN', 'R')
df_plot_c2 <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
pred_spatial_cell2))
names(df_plot_c2)[c(1,2,3)] <- c('EQID', 'WESN', 'R')
# # EQ 138
ref_id <- 138
df_plot_cb <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_cb2 <- df_plot_c2[df_plot_c2$EQID == ref_id
& df_plot_c2$WESN != 0,]
ggplot()+
geom_line(df_plot_cb, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_cb2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
ylim(-0.75, 3) +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
Here we plot the standard deviations associated with the predictions, which is a measure of predictive epistemic uncertainty.
# # EQ 138
ref_id <- 138
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_f2 <- df_plot_f[df_plot_f$EQID == ref_id
& df_plot_f$WESN != 0,]
p138 <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = sd)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_f2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN)), linetype = 'twodash') +
scale_x_continuous(trans='log10') +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
# EQ 139
ref_id <- 139
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_f2 <- df_plot_f[df_plot_f$EQID == ref_id
& df_plot_f$WESN != 0,]
p139 <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = sd)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_f2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN)), linetype = 'twodash') +
scale_x_continuous(trans='log10') +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
# EQ 140
ref_id <- 140
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_f2 <- df_plot_f[df_plot_f$EQID == ref_id
& df_plot_f$WESN != 0,]
p140 <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = sd)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_f2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN)), linetype = 'twodash') +
scale_x_continuous(trans='log10') +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
grid.arrange(p138, p139, p140, nrow = 2)
## plots for paper
wid <- 8
asp <- 0.8
df_plot <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, pred))
names(df_plot)[c(1,2,3,4)] <- c('EQID', 'WESN', 'R', 'Y_gwr')
df_plot_c <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, pred_spatial_cell))
names(df_plot_c)[c(1,2,3,4)] <- c('EQID', 'WESN', 'R', 'Y_gwr')
# # EQ 138
ref_id <- 138
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
p <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("red", "blue", "magenta", "cyan"),
labels=c("West", "East", "South", "North")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(log[10]," PGA")))) +
ylim(-0.75, 3) +
annotate(geom="label", x = 100, y = 3, label = "EQ 1", size = 10, fill = "white") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.2),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_eq1.pdf'), p,
width = wid, height = asp * wid)
# # EQ 139
ref_id <- 139
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
p <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("red", "blue", "magenta", "cyan"),
labels=c("West", "East", "South", "North")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(log[10]," PGA")))) +
ylim(-0.75, 3) +
annotate(geom="label", x = 100, y = 3, label = "EQ 2", size = 10, fill = "white") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.2),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_eq2.pdf'), p,
width = wid, height = asp * wid)
# # EQ 140
ref_id <- 140
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
p <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("red", "blue", "magenta", "cyan"),
labels=c("West", "East", "South", "North")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(log[10]," PGA")))) +
ylim(-0.75, 3) +
annotate(geom="label", x = 100, y = 3, label = "EQ 3", size = 10, fill = "white") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.2),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_eq3.pdf'), p,
width = wid, height = asp * wid)
tmp <- unique(co_eq_pred)
p <- ggplot() + theme_bw() + gg(mesh) +
geom_point(data = as.data.frame(co_stat), aes(x=longitude_st, y=latitude_st), fill= 'firebrick3',
size = 2, shape = 24, stroke = 0.7)+
geom_point(data = as.data.frame(co_eq), aes(x=longitude_ev, y=latitude_ev), fill= 'cadetblue2',
size = 2, shape = 21, stroke = 1.5)+
geom_point(data = as.data.frame(co_stat_pred), aes(x=X_stat,y=Y_stat), color = "cyan", size = 1) +
geom_point(data = as.data.frame(tmp), aes(x=X_ev,y=Y_ev), color = "red", size = 1) +
labs(x="X (km)", y="Y (km)") +
theme(axis.title = element_text(size=30), axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20)) +
annotate(geom="label", x = tmp$X_ev[1] + 70, y = tmp$Y_ev[1] + 50, label = "EQ 1", size = 4, fill = "white") +
annotate(geom="label", x = tmp$X_ev[2] + 70, y = tmp$Y_ev[2] + 50, label = "EQ 2", size = 4, fill = "white") +
annotate(geom="label", x = tmp$X_ev[3] + 70, y = tmp$Y_ev[3] + 50, label = "EQ 3", size = 4, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_mesh_pred.pdf'), p,
width = wid, height = asp * wid)10.1 Predictions 2
data_pred <- read.csv("/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/data_prediction_line.csv")
data_dm_comb <- rstan::read_rdump('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/dm_combined_line_25x25.Rdata')
dm_sparse_c <- as(data_dm_comb$RC,"dgCMatrix")
load(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY','RESULTS', 'results_pred_mswgr_line.Rdata'))
attach(data_pred)
b1p = (mag-mh)*(mag<=mh)
b2p = (mag-mh)*(mag>mh)
c1p = (mag-mref)*log10(sqrt(JB_complete^2+h^2))
c2p = log10(sqrt(JB_complete^2+h^2))
c3p = sqrt(JB_complete^2+h^2)
f1p = as.numeric(fm_type_code == "SS")
f2p = as.numeric(fm_type_code == "TF")
kp = log10(vs30/800)*(vs30<=1500)+log10(1500/800)*(vs30>1500)
yp = log10(rotD50_pga)
detach(data_pred)
data_pred_reg <- data.frame(Y = yp,
M1 = b1p,
M2 = b2p,
MlnR = c1p,
lnR = c2p,
R = c3p,
Fss = f1p,
Frv = f2p,
lnVS = kp,
#eq = data_pred$EQID,
eq = NA,
stat = data_pred$STATID,
intercept = 1,
resid = NA
)
data_pred_reg$idx_cell <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
data_pred_reg$idx <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
co_eq_pred <- data_pred[,c(13,14)]
co_stat_pred <- data_pred[,c(15,16)]
co_eq_stat_pred <- cbind(rowMeans(cbind(co_eq_pred$X_ev, co_stat_pred$X_stat)),
rowMeans(cbind(co_eq_pred$Y_ev, co_stat_pred$Y_stat)))
### st station index to NA, to avid using it in calculating uncertainty
data_pred_reg_ns <- data.frame(Y = yp,
M1 = b1p,
M2 = b2p,
MlnR = c1p,
lnR = c2p,
R = c3p,
Fss = f1p,
Frv = f2p,
lnVS = kp,
#eq = data_pred$EQID,
eq = NA,
#stat = data_pred$STATID,
stat = NA,
intercept = 1,
resid = NA
)
data_pred_reg_ns$idx_cell <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
data_pred_reg_ns$idx <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))A_stat_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_stat_pred))
A_eq_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_pred))
idx.stat_pred <- inla.spde.make.index("idx.stat",spde_stat$n.spde)
idx.eq_pred <- inla.spde.make.index("idx.eq",spde_eq$n.spde)
A_lnr_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_pred),
weights = data_pred_reg$lnR)
A_r_es_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_stat_pred),
weights = data_pred_reg$R)
idx.lnr_pred <- inla.spde.make.index("idx.lnr",spde_lnR$n.spde)
idx.r_pred <- inla.spde.make.index("idx.r",spde_R$n.spde)
A_vs_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_stat_pred),
weights = data_pred_reg$lnVS)
idx.vs_pred <- inla.spde.make.index("idx.vs",spde_vs$n.spde)
# create the stack
stk_spatial_full_pred <- inla.stack(
data = list(y = NA),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_eq_pred, A_lnr_pred, A_r_es_pred, A_stat_pred, A_vs_pred, 1),
effects = list(idx.eq = idx.eq_pred,
idx.lnr = idx.lnr_pred,
idx.r = idx.r_pred,
idx.stat = idx.stat_pred,
idx.vs = idx.vs_pred,
data_pred_reg
),
tag = 'model_spatial_full_pred')
stk_spatial_full_comb <- inla.stack(stk_spatial_full, stk_spatial_full_pred)
# create the stack
stk_spatial_full_pred_ns <- inla.stack(
data = list(y = NA),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_eq_pred, A_lnr_pred, A_r_es_pred, A_stat_pred, A_vs_pred, 1),
effects = list(idx.eq = idx.eq_pred,
idx.lnr = idx.lnr_pred,
idx.r = idx.r_pred,
idx.stat = idx.stat_pred,
idx.vs = idx.vs_pred,
data_pred_reg_ns
),
tag = 'model_spatial_full_pred_ns')
stk_spatial_full_comb_ns <- inla.stack(stk_spatial_full, stk_spatial_full_pred_ns)data_comb <- rbind(data_reg, data_pred_reg)
fit_inla_pred <- inla(form,
data = data_comb,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.predictor = list(compute = TRUE)
)
idx_pred <- inla.stack.index(stk_spatial_full_comb, 'model_spatial_full_pred')$data
pred <- fit_inla_pred$summary.fitted.values[idx_pred,]
data_comb <- rbind(data_reg, data_pred_reg_ns)
fit_inla_pred <- inla(form,
data = data_comb,
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.predictor = list(compute = TRUE)
)
pred_ns <- fit_inla_pred$summary.fitted.values[idx_pred,]
rm(fit_inla_pred)form_spatial_cell_c <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx_cell, model = "z", Z = dm_sparse_c, hyper = prior_prec_cell)
fit_inla_spatial_cell_pred <- inla(form_spatial_cell_c,
data = inla.stack.data(stk_spatial_full_comb),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)## Warning in inla.model.properties.generic(inla.trim.family(model), mm[names(mm) == : Model 'z' in section 'latent' is marked as 'experimental'; changes may appear at any time.
## Use this model with extra care!!! Further warnings are disabled.idx_pred <- inla.stack.index(stk_spatial_full_comb, 'model_spatial_full_pred')$data
pred_spatial_cell <- fit_inla_spatial_cell_pred$summary.fitted.values[idx_pred,]
### no stat
fit_inla_spatial_cell_pred <- inla(form_spatial_cell_c,
data = inla.stack.data(stk_spatial_full_comb_ns),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb_ns),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family = list(hyper = list(prec = prior_prec_phiSS_lg)),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)
pred_spatial_cell_ns <- fit_inla_spatial_cell_pred$summary.fitted.values[idx_pred,]
rm(fit_inla_spatial_cell_pred)df_plot <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
pred))
names(df_plot)[c(1,2,3)] <- c('EQID', 'WESN', 'R')
df_plot_ns <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
pred_ns))
names(df_plot_ns)[c(1,2,3)] <- c('EQID', 'WESN', 'R')
df_plot_c <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, sqrt(prediction$s0 + 1/fit_inla_resid$summary.hyperpar$mean[3]),
pred_spatial_cell))
names(df_plot_c)[c(1,2,3,4,5)] <- c('EQID', 'WESN', 'R', 'Y_gwr', 'SD_gwr')
df_plot_c_ns <- as.data.frame(cbind(data_pred$EQID, data_pred$WESN, data_pred$JB_complete,
prediction$y0, sqrt(prediction$s0),
pred_spatial_cell_ns))
names(df_plot_c_ns)[c(1,2,3,4,5)] <- c('EQID', 'WESN', 'R', 'Y_gwr', 'SD_gwr')
# # EQ 138
ref_id <- 138
df_plot2 <- df_plot[df_plot$EQID == ref_id
& df_plot$WESN == 1,]
df_plot2_ns <- df_plot_ns[df_plot_ns$EQID == ref_id
& df_plot_ns$WESN == 1,]
df_plot_c2 <- df_plot_c[df_plot_c$EQID == ref_id
& df_plot_c$WESN != 0,]
df_plot_c_ns2 <- df_plot_c_ns[df_plot_c_ns$EQID == ref_id
& df_plot_c_ns$WESN != 0,]
ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
ylim(-0.75, 3) +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = sd)) +
geom_line(df_plot2_ns, mapping = aes(x = R, y = sd), linetype = 'dotted') +
geom_line(df_plot_c2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c_ns2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN)), linetype = 'dotted') +
geom_line(df_plot_c2, mapping = aes(x = R, y = SD_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
labs(title = paste0("EQ ",ref_id)) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
tmp <- unique(co_eq_pred)
p <- ggplot() + theme_bw() + gg(mesh) +
geom_point(data = as.data.frame(co_stat), aes(x=longitude_st, y=latitude_st), fill= 'firebrick3',
size = 2, shape = 24, stroke = 0.7)+
geom_point(data = as.data.frame(co_eq), aes(x=longitude_ev, y=latitude_ev), fill= 'cadetblue2',
size = 2, shape = 21, stroke = 0.7)+
geom_line(data = as.data.frame(data_pred[data_pred$WESN == 1,c(15,16)]), aes(x=X_stat,y=Y_stat), color = "blue", size = 1.5) +
geom_line(data = as.data.frame(data_pred[data_pred$WESN == 2,c(15,16)]), aes(x=X_stat,y=Y_stat), color = "red", size = 1.5) +
geom_point(data = as.data.frame(tmp), aes(x=X_ev,y=Y_ev), fill= 'red',
size = 2, shape = 21, stroke = 1.5) +
labs(x="X (km)", y="Y (km)") +
theme(axis.title = element_text(size=30), axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20))
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_mesh_pred_line.pdf'), p,
width = wid, height = asp * wid)
p <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = mean)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = mean, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = Y_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("blue", "red"),
labels=c("North", "South")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(log[10]," PGA")))) +
ylim(0., 3) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.2),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_line.pdf'), p,
width = wid, height = asp * wid)## Warning: Removed 11 row(s) containing missing values (geom_path).p <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = sd)) +
geom_line(df_plot_c2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c2, mapping = aes(x = R, y = SD_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("blue", "red"),
labels=c("North", "South")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(psi[mu])))) +
ylim(0., 0.5) +
annotate(geom="text", x = 80, y = 0.45, label = expression(atop(paste("with ",phi[S2S]))),
size = 8) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.8),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_line_sd.pdf'), p,
width = wid, height = asp * wid)## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'p <- ggplot()+
geom_line(df_plot2, mapping = aes(x = R, y = sd)) +
geom_line(df_plot2_ns, mapping = aes(x = R, y = sd), linetype = 'dotted') +
geom_line(df_plot_c2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c_ns2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN)), linetype = 'dotted') +
geom_line(df_plot_c2, mapping = aes(x = R, y = SD_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("blue", "red"),
labels=c("North", "South")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(psi[mu])))) +
ylim(0., 0.5) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.8),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_line_sd_ns.pdf'), p,
width = wid, height = asp * wid)
p <- ggplot()+
geom_line(df_plot2_ns, mapping = aes(x = R, y = sd)) +
geom_line(df_plot_c_ns2, mapping = aes(x = R, y = sd, group = factor(WESN),
color = factor(WESN))) +
geom_line(df_plot_c_ns2, mapping = aes(x = R, y = SD_gwr, group = factor(WESN),
color = factor(WESN)), linetype = 'dashed') +
scale_x_continuous(trans='log10') +
annotation_logticks(sides = "b") +
scale_color_manual(values=c("blue", "red"),
labels=c("North", "South")
) +
guides(color=guide_legend(title="Direction")) +
xlab(expression(atop(paste(R[JB]," (km)")))) + ylab(expression(atop(paste(psi[mu])))) +
ylim(0., 0.5) +
annotate(geom="text", x = 80, y = 0.45, label = expression(atop(paste("without ",phi[S2S]))),
size = 8) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.1, 0.8),
legend.title = element_text(size = 20)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_pred_line_sd_ns2.pdf'), p,
width = wid, height = asp * wid)## Warning in is.na(x): is.na() applied to non-(list or vector) of type
## 'expression'10.2 Predictive Uncertainty
In the previous section, we had plotted the uncertainty associated with prediction. This prediction did not include \(\tau\), but all other terms (i.e. uncertainty due to fixed effects, site terms \(\delta S2S\), and spatial effects). Here we look at the full predictive uncertainty, which includes all terms. This is generally not something we need for PSHA, since the full predictive uncerainty is decomposed into aleatory vaiabiliy and epistemic uncertainty; the former is integrated out, the latter generally is not (and is used to assess a distribution of hazard curves; simply put, it is put on a logic tree). Looking a the full predictive distribution gives a good measure how an ergodic model compares with the nonergodic models.
In INLA, the posterior distribution associated with predictions is calculated when the value of the target variable is NA. This distribution represents the uncertainty associated with the mean prediction, and does not include the ``noise’’ term. One way to include this is to sample from the posterior distribution of the model, and then calculate summary statistics (like intervals, or standard deviations). This was done in the cross-validation section. Here, we use a trick taken from https://julianfaraway.github.io/brinla/examples/chicago.html (???). We fix the observation precision to a high value, and add a random effect corresponding to each observation. Since the observation standard deviation is now very small, the new random effect will absorb the noise term. Because it is a random effect, it is included in the prediction.
We predict for a new set of coordinates. These are not necessarily realistic event or station locations, but are chosen so that they are far away from observed data.
data_pred <- read.csv("/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/data_prediction_pt.csv")
data_dm_comb <- rstan::read_rdump('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/dm_combined_pt_25x25.Rdata')
dm_sparse_c <- as(data_dm_comb$RC,"dgCMatrix")
attach(data_pred)
b1p = (mag-mh)*(mag<=mh)
b2p = (mag-mh)*(mag>mh)
c1p = (mag-mref)*log10(sqrt(JB_complete^2+h^2))
c2p = log10(sqrt(JB_complete^2+h^2))
c3p = sqrt(JB_complete^2+h^2)
f1p = as.numeric(fm_type_code == "SS")
f2p = as.numeric(fm_type_code == "TF")
kp = log10(vs30/800)*(vs30<=1500)+log10(1500/800)*(vs30>1500)
yp = log10(rotD50_pga)
detach(data_pred)
data_pred_reg <- data.frame(Y = yp,
M1 = b1p,
M2 = b2p,
MlnR = c1p,
lnR = c2p,
R = c3p,
Fss = f1p,
Frv = f2p,
lnVS = kp,
eq = data_pred$EQID,
stat = data_pred$STATID,
intercept = 1
)
data_pred_reg$idx_cell <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
data_pred_reg$resid <- NA
# add index corresponding to noise term for each data and prediction point
data_pred_reg$idx <- (n_rec + 1):(n_rec + length(data_pred_reg$Y))
co_eq_pred <- data_pred[,c(13,14)]
co_stat_pred <- data_pred[,c(15,16)]
co_eq_stat_pred <- cbind(rowMeans(cbind(co_eq_pred$X_ev, co_stat_pred$X_stat)),
rowMeans(cbind(co_eq_pred$Y_ev, co_stat_pred$Y_stat)))
ggplot() + theme_bw() + gg(mesh) +
geom_point(data = as.data.frame(co_stat_pred), aes(x=X_stat,y=Y_stat), color = set1[1]) +
geom_point(data = as.data.frame(co_eq_pred), aes(x=X_ev,y=Y_ev), color = set1[2]) +
labs(x="X (km)", y="Y (km)") +
theme(axis.title = element_text(size=30), axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20))
Calculate projection matrices, and prediction stack.
A_stat_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_stat_pred))
A_eq_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_pred))
idx.stat_pred <- inla.spde.make.index("idx.stat",spde_stat$n.spde)
idx.eq_pred <- inla.spde.make.index("idx.eq",spde_eq$n.spde)
A_lnr_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_pred),
weights = data_pred_reg$lnR)
A_r_es_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_eq_stat_pred),
weights = data_pred_reg$R)
idx.lnr_pred <- inla.spde.make.index("idx.lnr",spde_lnR$n.spde)
idx.r_pred <- inla.spde.make.index("idx.r",spde_R$n.spde)
A_vs_pred <- inla.spde.make.A(mesh, loc = as.matrix(co_stat_pred),
weights = data_pred_reg$lnVS)
idx.vs_pred <- inla.spde.make.index("idx.vs",spde_vs$n.spde)
# create the stack
stk_spatial_full_pred <- inla.stack(
data = list(y = NA),
#A = list(A_eq, A_eq, A_eq, A_stat, A_stat, 1),
A = list(A_eq_pred, A_lnr_pred, A_r_es_pred, A_stat_pred, A_vs_pred, 1),
effects = list(idx.eq = idx.eq_pred,
idx.lnr = idx.lnr_pred,
idx.r = idx.r_pred,
idx.stat = idx.stat_pred,
idx.vs = idx.vs_pred,
data_pred_reg
),
tag = 'model_spatial_full_pred')
stk_spatial_full_comb <- inla.stack(stk_spatial_full, stk_spatial_full_pred)
idx_pred <- inla.stack.index(stk_spatial_full_comb, 'model_spatial_full_pred')$data
data_comb <- rbind(data_reg, data_pred_reg)
# create stack withBelow, we refit the models with the fixed likelihood and the new random effect to get the full predictive uncertainty distribution. We also fit a model where the cells are correlated. This model is based on an autoregressive model (Ver Hoef et al. 2018), where the correlation is determined based on neighborhood structure. For a fixed value of the correlation coefficient \(\rho\), under a conditional autoregressive (CAR) model the precision matrix becomes
\[ \Sigma^{-1} = \textbf{I} - \rho \textbf{W} \] where \(\textbf{I}\) is the identity matrix, and \(\textbf{W}\) is the adjacency matrix describing the neighborhood structure. In this case, if two cells share a border, the entry in the adjacency matix is 1, otherwise it is zero (rook adjacency).
# ergodic base model
form2 <- Y ~ M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx, model = "iid",hyper = prior_prec_phiSS_lg)
fit_inla_pred2 <- inla(form2,
data = data_comb,
family="gaussian",
control.fixed = prior.fixed,
control.family=list(initial=12,fixed=TRUE),
control.predictor = list(compute = TRUE)
)
pred2 <- fit_inla_pred2$summary.fitted.values[idx_pred,]
# spatial model with event and station constants and cells
form_spatial_cell_c2 <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx_cell, model = "z", Z = dm_sparse_c, hyper = prior_prec_cell) +
f(idx, model = "iid", hyper = prior_prec_phiSS_lg)
fit_inla_spatial_cell_pred2 <- inla(form_spatial_cell_c2,
data = inla.stack.data(stk_spatial_full_comb),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family=list(initial=12,fixed=TRUE),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)
pred_spatial_cell2 <- fit_inla_spatial_cell_pred2$summary.fitted.values[idx_pred,]
# spatial model with event and station constants
form_spatial2 <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx, model = "iid",hyper = prior_prec_phiSS_lg)
fit_inla_spatial_pred2 <- inla(form_spatial2,
data = inla.stack.data(stk_spatial_full_comb),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family=list(initial=12,fixed=TRUE),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)
pred_spatial2 <- fit_inla_spatial_pred2$summary.fitted.values[idx_pred,]
##### correlation of cells
adj <- read.csv("/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY/DATA/cells_adj.csv")
cell_idx <- data_dm_comb$cell_idx
rho <- 0.2
tmp <- diag(length(adj[,1])) - rho * adj
cmat <- as(as.matrix(tmp[cell_idx, cell_idx]),"dgCMatrix")
form_spatial_cell_corr_c2 <- y ~ 0 + intercept + M1 + M2 + lnR + MlnR + R + Fss + Frv + lnVS +
f(eq, model = "iid", hyper = prior_prec_tau_lg) +
f(stat, model = "iid",hyper = prior_prec_phiS2S_lg) +
f(idx.stat, model = spde_stat) + f(idx.eq, model = spde_eq) +
f(idx_cell, model = "z", Z = dm_sparse_c, hyper = prior_prec_cell, Cmatrix = cmat) +
f(idx, model = "iid", hyper = prior_prec_phiSS_lg)
fit_inla_spatial_cell_corr_pred2 <- inla(form_spatial_cell_corr_c2,
data = inla.stack.data(stk_spatial_full_comb),
control.predictor = list(A = inla.stack.A(stk_spatial_full_comb),
compute = TRUE),
family="gaussian",
control.fixed = prior.fixed,
control.family=list(initial=12,fixed=TRUE),
control.inla = list(int.strategy = "eb", strategy = "gaussian")
)
pred_spatial_cell_corr2 <- fit_inla_spatial_cell_corr_pred2$summary.fitted.values[idx_pred,]df_plot <- data_pred_reg
df_plot$mean_inla2 <- pred2$mean
df_plot$mean_spatial_cell2 <- pred_spatial_cell2$mean
df_plot$mean_spatial_cell_corr2 <- pred_spatial_cell_corr2$mean
df_plot$mean_spatial2 <- pred_spatial2$mean
df_plot$sd_inla2 <- pred2$sd
df_plot$sd_spatial_cell2 <- pred_spatial_cell2$sd
df_plot$sd_spatial_cell_corr2 <- pred_spatial_cell_corr2$sd
df_plot$sd_spatial2 <- pred_spatial2$sd
ggplot(df_plot) +
geom_line(aes(x = R, y = mean_inla2)) +
geom_line(aes(x = R, y = mean_spatial2), color = 'blue') +
geom_line(aes(x = R, y = mean_spatial_cell2), color = 'red') +
geom_line(aes(x = R, y = mean_spatial_cell_corr2), color = 'darkred', linetype = 'dashed') +
scale_x_continuous(trans='log10') +
labs(title = "Prediction",
subtitle = "",
x = "R(km)", y = "mean") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
ggplot(df_plot) +
geom_line(aes(x = R, y = sd_inla2)) +
geom_line(aes(x = R, y = sd_spatial2), color = 'blue') +
geom_line(aes(x = R, y = sd_spatial_cell2), color = 'red') +
geom_line(aes(x = R, y = sd_spatial_cell_corr2), color = 'darkred', linetype = 'dashed') +
scale_x_continuous(trans='log10') +
labs(title = "Prediction",
subtitle = "Full predictive uncertainty",
x = "R(km)", y = "sd") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20)
)
11 Some Remarks
According to Caramenti et al. (2020) and Lanzano et al. (2021), MS-GWR has some advantages
- It is easier to understand
- maybe true; a local weighted model is easy to understand, but a GP (basically a large multi-normal is not hard)
- it is non-parametric (meaning no distribution for the response is assumed)
- that is true, but kind of irrelevant, since for a GMM one must make some distributional assumption (a GMM provides the gistribution of GM, not just a median prediction)
Based on this investigation, GP based VCMs fit with INLA fit the data better (have lower generalization error and higher loglikelihood based on 5fold CV). Another bgh advantage is time; CV for MS-GWR takes mre time than for the 6 INLA models (four of which are spatial models).
12 Some additional plots
wid <- 8
asp <- 0.8
### phi SS
res <- rbind(1/sqrt(fit_inla_noranef$summary.hyperpar[,c(4,5,3)]),
1/sqrt(fit_inla$summary.hyperpar[1,c(4,5,3)]),
1/sqrt(fit_inla_spatial$summary.hyperpar[1,c(4,5,3)]),
1/sqrt(fit_inla_spatial_full$summary.hyperpar[1,c(4,5,3)]),
1/sqrt(fit_inla_spatial_cell$summary.hyperpar[1,c(4,5,3)]),
1/sqrt(fit_inla_spatial_full_b$summary.hyperpar[1,c(4,5,3)]),
1/sqrt(fit_inla_resid$summary.hyperpar[1,c(4,5,3)])
)
colnames(res) = c("E", "L", "U")
res$model <- c("M1", "M3", "M4", "M7", "M9", "M10", "MS-WGR B")
res$model <- factor(res$model, levels = res$model)
p <- ggplot(res, aes(x = model, y = E)) +
geom_point(size = 3) +
geom_pointrange(aes(ymax = U, ymin = L)) +
ylim(0,0.4) +
xlab(NULL) + ylab(expression(atop(paste(phi[SS])))) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_phiss.pdf'), p,
width = wid, height = asp * wid)
### tau
res <- rbind(c(NA, NA, NA),
1/sqrt(fit_inla$summary.hyperpar[2,c(4,5,3)]),
1/sqrt(fit_inla_spatial$summary.hyperpar[2,c(4,5,3)]),
1/sqrt(fit_inla_spatial_full$summary.hyperpar[2,c(4,5,3)]),
1/sqrt(fit_inla_spatial_cell$summary.hyperpar[2,c(4,5,3)]),
c(NA, NA, NA),
1/sqrt(fit_inla_resid$summary.hyperpar[2,c(4,5,3)])
)
colnames(res) = c("E", "L", "U")
res$model <- c("M1", "M3", "M4", "M7", "M9", "M10", "MS-WGR B")
res$model <- factor(res$model, levels = res$model)
p <- ggplot(res, aes(x = model, y = E)) +
geom_point(size = 3) +
geom_pointrange(aes(ymax = U, ymin = L)) +
xlab(NULL) + ylab(expression(atop(paste(tau)))) +
ylim(0,0.25) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_tau.pdf'), p,
width = wid, height = asp * wid)## Warning: Removed 2 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_pointrange).### phi s2s
res <- rbind(c(NA, NA, NA),
1/sqrt(fit_inla$summary.hyperpar[3,c(4,5,3)]),
1/sqrt(fit_inla_spatial$summary.hyperpar[3,c(4,5,3)]),
1/sqrt(fit_inla_spatial_full$summary.hyperpar[3,c(4,5,3)]),
1/sqrt(fit_inla_spatial_cell$summary.hyperpar[3,c(4,5,3)]),
c(NA, NA, NA),
1/sqrt(fit_inla_resid$summary.hyperpar[3,c(4,5,3)])
)
colnames(res) = c("E", "L", "U")
res$model <- c("M1", "M3", "M4", "M7", "M9", "M10", "MS-WGR B")
res$model <- factor(res$model, levels = res$model)
p <- ggplot(res, aes(x = model, y = E)) +
geom_point(size = 3) +
geom_pointrange(aes(ymax = U, ymin = L)) +
xlab(NULL) + ylab(expression(atop(paste(phi[S2S])))) +
ylim(0,0.25) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_phis2s.pdf'), p,
width = wid, height = asp * wid)## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_pointrange).# # all standard deviations
sd_list <- c("sigma", "phi_ss", "tau", "phi_s2s")
res <- rbind(1/sqrt(fit_inla_noranef$summary.hyperpar[,c(4,5,3)]),c(NA,NA,NA),c(NA,NA,NA),c(NA,NA,NA))
res <- rbind(res, c(NA,NA,NA), 1/sqrt(fit_inla$summary.hyperpar[,c(4,5,3)]))
res <- rbind(res, c(NA,NA,NA), 1/sqrt(fit_inla_spatial$summary.hyperpar[c(1,2,3),c(4,5,3)]))
res <- rbind(res, c(NA,NA,NA), 1/sqrt(fit_inla_spatial_full$summary.hyperpar[c(1,2,3),c(4,5,3)]))
res <- rbind(res, c(NA,NA,NA), 1/sqrt(fit_inla_spatial_cell$summary.hyperpar[c(1,2,3),c(4,5,3)]))
res <- rbind(res, 1/sqrt(fit_inla_spatial_full_b$summary.hyperpar[1,c(4,5,3)]),c(NA,NA,NA),c(NA,NA,NA),c(NA,NA,NA))
res <- rbind(res, c(0.3008844, NA, NA), c(NA,NA,NA),c(NA,NA,NA),c(NA,NA,NA))
res <- rbind(res, c(NA,NA,NA),1/sqrt(fit_inla_resid$summary.hyperpar[c(1,2,3),c(4,5,3)]))
colnames(res) = c("E", "L", "U")
n.sd <- length(sd_list)
model_list <- c("Model 1", "Model 3", "Model 4", "Model 7", "Model 9", "Model 10",
"MS-GWR", "MS-GWR B")
res$model = factor(rep(model_list, each=n.sd),
levels = model_list)
res$sd = factor(rep(sd_list, length(model_list)), levels = sd_list)
p <- ggplot(res, aes(x = model, y = E, color = sd)) +
geom_point(size = 3) +
geom_pointrange(aes(ymax = U, ymin = L)) +
xlab(NULL) + ylab("posterior") +
scale_color_manual(values=c("black", "red", "blue", "magenta"),
labels=c(expression(atop(paste(sigma))),
expression(atop(paste(phi[SS]))),
expression(atop(paste(tau))),
expression(atop(paste(phi[S2S]))))
) +
ylim(0,0.4) +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
legend.text = element_text(size = 20),
legend.position = c(0.5, 0.8),
legend.box.background = element_rect(colour = "black"),
axis.text.x = element_text(angle = 45, hjust=1)
) +
guides(color = guide_legend(nrow = 2, byrow = TRUE,title=NULL))
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_sd_all.pdf'), p,
width = wid, height = asp * wid)## Warning: Removed 14 rows containing missing values (geom_point).## Warning: Removed 14 rows containing missing values (geom_pointrange).## Warning: Removed 1 rows containing missing values (geom_segment).###########
# fixed effects plots
res <- fit_inla$summary.fixed[,c(4,3,5)]
res <- rbind(res, fit_inla_spatial$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial_full$summary.fixed[,c(4,3,5)])
res <- rbind(res, fit_inla_spatial_cell$summary.fixed[,c(4,3,5)])
colnames(res) = c("E", "L", "U")
rownames(res)=NULL
n.covar = nrow(fit_inla$summary.fixed)
model_list <- c("M3", "M4", "M7", "M9")
res$model = factor(rep(model_list, each=n.covar),
levels = model_list)
covar_list <- rownames(fit_inla$summary.fixed)
covar_list <- c("a", "b[1]", "b[2]", "c[2]", "c[1]", "c[3]", "f[1]", "f[2]", "k")
res$covar = factor(rep(covar_list, length(model_list)),
levels = covar_list)
p <- ggplot(res, aes(x = model, y = E)) +
facet_wrap(~covar, scales = "free_y", labeller = label_parsed) +
geom_point(size = 3) +
geom_pointrange(aes(ymax = U, ymin = L)) +
xlab(NULL) + ylab(NULL) +
theme(
axis.title = element_text(size = 18),
axis.text = element_text(size = 12),
plot.title = element_text(size = 30),
strip.text.x = element_text(size = 12)
)
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_fixed_effects.pdf'), p,
width = wid, height = asp * wid)#### full spatal model
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.lnr$mean +
fit_inla_spatial_full$summary.fixed$mean[4])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[2])))
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = c(0.1, 0.2)
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 7", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_gs_full.pdf'), p,
width = wid, height = asp * wid)
#### full spatial model no random effects
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.lnr$mean +
fit_inla_spatial_full_b$summary.fixed$mean[4])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[2])))
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 10", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_gs_full_b.pdf'), p,
width = wid, height = asp * wid)
#### full spatial model no random effects
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_cell$summary.random$idx.lnr$mean +
fit_inla_spatial_full_cell$summary.fixed$mean[4])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[2])))
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 8", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_gs_full_cell.pdf'), p,
width = wid, height = asp * wid)
#### ms-wgr
df_plot <- as.data.frame(cbind(co_df_utm, beta_c2))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=beta_c2))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[2])))
,limits = c(-1.8,-1.), breaks = c(-1.8, -1.4, -1)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "MS-GWR", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_gs_msgwr.pdf'), p,
width = wid, height = asp * wid)#### full spatal model
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.vs$mean +
fit_inla_spatial_full$summary.fixed$mean[9])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "k"
,limits = c(-1.,0.00), breaks = c(-1., -0.5, 0.)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = c(0.1, 0.2)
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 7", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_vs_full.pdf'), p,
width = wid, height = asp * wid)
#### full spatial model no random effects
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.vs$mean +
fit_inla_spatial_full_b$summary.fixed$mean[9])
field.proj[field.proj < -1] <- -1
field.proj[field.proj > 0] <- 0
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "k"
,limits = c(-1.,0.00), breaks = c(-1., -0.5, 0.)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 10", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_vs_full_b.pdf'), p,
width = wid, height = asp * wid)
#### full spatial model no random effects
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_cell$summary.random$idx.vs$mean +
fit_inla_spatial_full_cell$summary.fixed$mean[9])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "k"
,limits = c(-1.,0.00), breaks = c(-1., -0.5, 0.)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 8", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_vs_full_cell.pdf'), p,
width = wid, height = asp * wid)
#### ms-gwr
df_plot <- as.data.frame(cbind(co_df_utm, beta_k))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=beta_k))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = "k"
,limits = c(-1., 0.00), breaks = c(-1., -0.5, 0.)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "MS-GWR", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_vs_msgwr.pdf'), p,
width = wid, height = asp * wid)#### full spatial model
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full$summary.random$idx.r$mean +
fit_inla_spatial_full$summary.fixed$mean[6])
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
field.proj[field.proj < -0.01] <- -0.01
field.proj[field.proj > 0.002] <- 0.002
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[3])))
,limits = c(-0.01,0.002), breaks = c(-0.01, -0.004, 0.002)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = c(0.13, 0.23)
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 7", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_linr_full.pdf'), p,
width = wid, height = asp * wid)
#### full spatial model no random effects
field.proj <- inla.mesh.project(proj, fit_inla_spatial_full_b$summary.random$idx.r$mean +
fit_inla_spatial_full_b$summary.fixed$mean[6])
field.proj[field.proj < -0.01] <- -0.01
field.proj[field.proj > 0.002] <- 0.002
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[3])))
,limits = c(-0.01,0.002), breaks = c(-0.01, -0.004, 0.002)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 10", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_linr_full_b.pdf'), p,
width = wid, height = asp * wid)
#### full spatial model cells
ca <- fit_inla_spatial_full_cell$summary.random$idx_cell$mean[(n_rec + 1):(n_rec + n_cell)]
field.proj <- ca + fit_inla_spatial_full_cell$summary.fixed$mean[6]
field.proj[field.proj < -0.01] <- -0.01
field.proj[field.proj > 0.002] <- 0.002
df_plot <- data.frame(x = cells2$X0, y = cells2$X1,
z = field.proj)
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x, y=y, fill=z))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[3])))
,limits = c(-0.01,0.002), breaks = c(-0.01, -0.004, 0.002)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "Model 8", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_linr_full_cell.pdf'), p,
width = wid, height = asp * wid)
#### ms-gwr
field.proj <- beta_c3
field.proj[field.proj < -0.01] <- -0.01
field.proj[field.proj > 0.002] <- 0.002
df_plot <- as.data.frame(cbind(co_df_utm, field.proj))
p <- ggplot() +
geom_tile(df_plot, mapping = aes(x=x1, y=x2, fill=field.proj))+
scale_fill_gradientn(colours = c("darkblue", "dodgerblue1", "cadetblue2", "white"), name = expression(atop(paste(c[3])))
,limits = c(-0.01,0.002), breaks = c(-0.01, -0.004, 0.002)
) +
labs(x = "UTM X (KM)", y = "UTM Y (KM)") +
theme(
axis.title = element_text(size = 20),
axis.text = element_text(size = 14),
plot.title = element_text(size = 30),
plot.subtitle = element_text(size = 20),
legend.text = element_text(size = 20),
legend.title = element_text(size = 20),
legend.position = "none"
) +
annotate(geom="label", x = 600, y = 5200, label = "MS-GWR", size = 10, fill = "white")
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_spatial_linr_msgwr.pdf'), p,
width = wid, height = asp * wid)p <- ggplot() + theme_bw() + gg(mesh) +
geom_point(data = as.data.frame(co_stat), aes(x=longitude_st, y=latitude_st), fill= 'firebrick3',
size = 2, shape = 24, stroke = 0.7)+
geom_point(data = as.data.frame(co_eq), aes(x=longitude_ev, y=latitude_ev), fill= 'cadetblue2',
size = 2, shape = 21, stroke = 1.5)+
labs(x="X (km)", y="Y (km)") +
theme(axis.title = element_text(size=30), axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20))
ggsave(file.path('/Users/nico/GROUNDMOTION/PROJECTS/NONERGODIC/ITALY',
'RESULTS','PLOTS','plot_mesh.pdf'), p,
width = wid, height = asp * wid)ca <- fit_inla_spatial_full_cell$summary.random$idx_cell$mean[(n_rec + 1):(n_rec + n_cell)]
field.proj <- ca + fit_inla_spatial_full_cell$summary.fixed$mean[6]
hist(field.proj)
sum(field.proj > 0)## [1] 15## NULL
13 Acknowledgments
This tutorial was written using the package markdown (Allaire et al. 2019) and rendered with knitr (Xie 2014, 2015, 2021). The R Markdown Cookbook (Xie, Dervieux, and Riederer 2020) was a helpful resource for R-Markdown. The package inlabru (Bachl et al. 2019) was used to make the mesh plot.
References
Abrahamson, Norman, Nicolas Kuehn, Melanie Walling, and Niels Landwehr. 2019. “Probabilistic Seismic Hazard Analysis in California Using Nonergodic Ground Motion Models.” Bulletin of the Seismological Society of America 109 (4): 1235–49. https://doi.org/10.1785/0120190030.
Akaike, Hirotogu. 1973. “Information Theory and an Extension of the Maximum Likelihood Principle.” Second International Symposium on Information Theory, 267–81. https://doi.org/10.1007/978-1-4612-1694-0_15.
Allaire, J J, Jeffrey Horner, Yihui Xie, Vicent Marti, and Natacha Porte. 2019. markdown: Render Markdown with the C Library ’Sundown’. https://cran.r-project.org/package=markdown.
Anderson, J G, and J N Brune. 1999. “Probabilistic Seismic Hazard Analysis without the Ergodic Assumption.” Seismological Research Letters 70 (1): 19–28. https://doi.org/10.1785/gssrl.70.1.19.
Bachl, Fabian E., Finn Lindgren, David L. Borchers, and Janine B. Illian. 2019. “inlabru: an R package for Bayesian spatial modelling from ecological survey data.” Methods in Ecology and Evolution 10 (6): 760–66. https://doi.org/10.1111/2041-210X.13168.
Bakka, Haakon, Håvard Rue, Geir-Arne Fuglstad, Andrea Riebler, David Bolin, Janine Illian, Elias Krainski, Daniel Simpson, and Finn Lindgren. 2018. “Spatial modeling with R-INLA: A review.” Wiley Interdisciplinary Reviews: Computational Statistics 10 (6): e1443. https://doi.org/10.1002/wics.1443.
Bussas, Matthias, Christoph Sawade, Nicolas Kühn, Tobias Scheffer, and Niels Landwehr. 2017. “Varying-coefficient models for geospatial transfer learning.” Machine Learning 106 (9-10): 1419–40. https://doi.org/10.1007/s10994-017-5639-3.
Caramenti, L., A. Menafoglio, S. Sgobba, and G. Lanzano. 2020. “Multi-Source Geographically Weighted Regression for Regionalized Ground-Motion Models.” https://mox.polimi.it/publication-results/?id=917%5C&tipo=add_qmox.
Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “Stan : A Probabilistic Programming Language.” Journal of Statistical Software 76 (1): 1–32. https://doi.org/10.18637/jss.v076.i01.
Dawood, Haitham M, and Adrian Rodriguez-Marek. 2013. “A Method for Including Path Effects in Ground-Motion Prediction Equations: An Example Using the Mw 9.0 Tohoku Earthquake Aftershocks.” Bulletin of the Seismological Society of America 103 (2B): 1360–72. https://doi.org/10.1785/0120120125.
Fotheringham, A. S., C. Brunsdon, and M. Charlton. 2002. Geographically Weighted Regression: The Analysis of Spatially Varying Relationships. Wiley.
Franco-Villoria, Maria, Massimo Ventrucci, and Håvard Rue. 2018. “Bayesian varying coefficient models using PC priors,” 1–21. http://arxiv.org/abs/1806.02084.
Fuglstad, Geir-Arne, Daniel Simpson, Finn Lindgren, and Håvard Rue. 2019. “Constructing Priors that Penalize the Complexity of Gaussian Random Fields.” Journal of the American Statistical Association 114 (525): 445–52. https://doi.org/10.1080/01621459.2017.1415907.
Gelfand, Alan E, Hyon-Jung Kim, C. F Sirmans, and Sudipto Banerjee. 2003. “Spatial Modeling with Spatially Varying Coefficient Processes.” Journal of the American Statistical Association 98 (462): 387–96. https://doi.org/10.1198/016214503000170.
Krainski, Elias, Virgilio Gómez-Rubio, Haakon Bakka, Amanda Lenzi, Daniela Castro-Camilo, Daniel Simpson, Finn K. Lindgren, and Håvard Rue. 2019. Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA. Boca-Raton, FL: Chapman; Hall/CRC.
Kuehn, Nicolas. 2021. “Comparison of Bayesian Varying Coefficient Models for the Development of Nonergodic Ground-Motion Models.” Engrxiv, 1–25. https://doi.org/10.31224/osf.io/tjxa3.
Kuehn, Nicolas Martin, Norman A Abrahamson, and Melanie Anne Walling. 2019. “Incorporating Nonergodic Path Effects into the NGA-West2 Ground-Motion Prediction Equations.” Bulletin of the Seismological Society of America 109 (2): 575–85. https://doi.org/10.1785/0120180260.
Landwehr, Niels, Nicolas M. Kuehn, Tobias Scheffer, and Norman Abrahamson. 2016. “A Nonergodic Ground‐Motion Model for California with Spatially Varying Coefficients.” Bulletin of the Seismological Society of America 106 (6): 2574–83. https://doi.org/10.1785/0120160118.
Lanzano, Giovanni, Lucia Luzi, Francesca Pacor, Chiara Felicetta, Rodolfo Puglia, Sara Sgobba, and Maria D’Amico. 2019. “A Revised Ground‐Motion Prediction Model for Shallow Crustal Earthquakes in Italy.” Bulletin of the Seismological Society of America 109 (2): 525–40. https://doi.org/10.1785/0120180210.
Lanzano, Giovanni, Sara Sgobba, Luca Caramenti, and Alessandra Menafoglio. 2021. “Ground-Motion Model for Crustal Events in Italy by Applying the Multisource Geographically Weighted Regression (MS-GWR) Method.” Bulletin of the Seismological Society of America, August, 1–17. https://doi.org/10.1785/0120210044.
Lavrentiadis, Grigorios, Norman A. Abrahamson, and Nicolas M. Kuehn. 2021. “A Non-ergodic Effective Amplitude Ground-Motion Model for California,” June. http://arxiv.org/abs/2106.07834.
Lindgren, Finn, and Håvard Rue. 2015. “Bayesian Spatial Modelling with R - INLA.” Journal of Statistical Software 63 (19): 1–25. https://doi.org/10.18637/jss.v063.i19.
Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An explicit link between gaussian fields and gaussian markov random fields: The stochastic partial differential equation approach.” Journal of the Royal Statistical Society. Series B: Statistical Methodology 73 (4): 423–98. https://doi.org/10.1111/j.1467-9868.2011.00777.x.
Rasmussen, Carl Edward, and Hannes Nickisch. 2010. “Gaussian Processes for Machine Learning (GPML) Toolbox.” Journal of Machine Learning Research 11: 3011–5. http://eprints.pascal-network.org/archive/00007469/.
Rasmussen, C E, and C K I Williams. 2006. Gaussian Processes for Machine Learning. Cambridge: MIT Press. http://www.gaussianprocess.org/gpml/.
Rue, Håvard, Sara Martino, and Nicolas Chopin. 2009. “Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 71 (2): 319–92. https://doi.org/10.1111/j.1467-9868.2008.00700.x.
Rue, Håvard, Andrea Riebler, Sigrunn H. Sørbye, Janine B. Illian, Daniel P. Simpson, and Finn K. Lindgren. 2017. “Bayesian Computing with INLA: A Review.” Annual Review of Statistics and Its Application 4 (1): 395–421. https://doi.org/10.1146/annurev-statistics-060116-054045.
Simpson, Daniel, Håvard Rue, Andrea Riebler, Thiago G. Martins, and Sigrunn H. Sørbye. 2017. “Penalising Model Component Complexity: A Principled, Practical Approach to Constructing Priors.” Statistical Science 32 (1): 1–28. https://doi.org/10.1214/16-STS576.
Snelson, Edward, and Zoubin Ghahramani. 2005. “Sparse Gaussian Processes using Pseudo-inputs.” In Advances in Neural Information Processing Systems, edited by Y. Weiss, B. Schoelkopf, and J. Platt. MIT Press. http://papers.nips.cc/paper/2857-sparse-gaussian-processes-using-pseudo-inputs.pdf%0Ahttp://www.gatsby.ucl.ac.uk/$\sim$snelson/SPGP_up.pdf.
Spiegelhalter, D J, N G Best, B P Carlin, and A Van Der Linde. 2002. “Bayesian measures of model complexity and fit.” Journal of the Royal Statistical Society. Series B, Statistical Methodology, 583–639. http://www.jstor.org/stable/3088806.
Ver Hoef, Jay M., Erin E. Peterson, Mevin B. Hooten, Ephraim M. Hanks, and Marie-Josèe Fortin. 2018. “Spatial autoregressive models for statistical inference from ecological data.” Ecological Monographs 88 (1): 36–59. https://doi.org/10.1002/ecm.1283.
Villani, Manuela, and Norman A Abrahamson. 2015. “Repeatable Site and Path Effects on the Ground‐Motion Sigma Based on Empirical Data from Southern California and Simulated Waveforms from the CyberShake Platform.” Bulletin of the Seismological Society of America 105 (5): 2681–95. https://doi.org/10.1785/0120140359.
Walling, M, and N A Abrahamson. 2012. “Non-Ergodic Probabilistic Seismic Hazard Analyses.” 15th World Conference on Earthquake Engineering (15WCEE).
Watanabe, Sumio. 2013. “A widely applicable bayesian information criterion.” Journal of Machine Learning Research 14 (1): 867–97. http://arxiv.org/abs/1208.6338.
Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D Peng. Chapman; Hall/CRC. http://www.crcpress.com/product/isbn/9781466561595.
———. 2015. Dynamic Documents with R and knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.
———. 2021. knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.
Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. 1st Editio. Chapman; Hall/CRC.
Unversity of California, Los Angeles, kuehn@ucla.edu↩︎